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

emdata_transform.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 "emdata.h"
00037 #include "emfft.h"
00038 
00039 #include <cstring>
00040 #include <cstdio>
00041 
00042 #include  "gsl_sf_result.h"
00043 #include  "gsl_sf_bessel.h"
00044 #include <iostream>
00045 #include <algorithm>
00046 #include <vector>
00047 #include <utility>
00048 #include <cmath>
00049 #include "util.h"
00050 
00051 #ifdef EMAN2_USING_CUDA
00052 #include "cuda/cuda_processor.h"
00053 #endif
00054 
00055 using namespace EMAN;
00056 using namespace std;
00057 typedef vector< pair<float,int> > vp;
00058 
00059 #ifdef EMAN2_USING_CUDA
00060 
00061 #include "cuda/cuda_emfft.h"
00062 
00063 EMData *EMData::do_fft_cuda() const
00064 {
00065         ENTERFUNC;
00066 
00067         if ( is_complex() ) {
00068                 LOGERR("real image expected. Input image is complex image.");
00069                 throw ImageFormatException("real image expected. Input image is complex image.");
00070         }
00071 //      int nxreal = nx;
00072         int offset;
00073         int ndim = get_ndim();
00074         EMData* dat = new EMData();
00075 
00076         offset = 2 - nx%2;
00077         dat->set_size_cuda(nx+offset,ny, nz);
00078         float *d = dat->get_cuda_data();
00079         if ( ndim == 1 ) {
00080                 cuda_dd_fft_real_to_complex_nd(get_cuda_data(), d, nx, 1,1);
00081         } else if (ndim == 2) {
00082                 cuda_dd_fft_real_to_complex_nd(get_cuda_data(), d, ny, nx, 1);
00083         } else if (ndim == 3) {
00084                 cuda_dd_fft_real_to_complex_nd(get_cuda_data(), d, nz, ny, nx);
00085         } else throw ImageDimensionException("No cuda FFT support of images with dimensions exceeding 3");
00086 
00087         if (offset == 1) dat->set_fftodd(true);
00088         else             dat->set_fftodd(false);
00089 
00090         dat->set_fftpad(true);
00091         dat->set_complex(true);
00092         if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(true);
00093         dat->set_ri(true);
00094         dat->gpu_update();
00095 
00096         EXITFUNC;
00097         return dat;
00098 }
00099 
00100 EMData *EMData::do_ift_cuda(bool preserve_input) const
00101 {
00102         ENTERFUNC;
00103 
00104         if (!is_complex()) {
00105                 LOGERR("complex image expected. Input image is real image.");
00106                 throw ImageFormatException("complex image expected. Input image is real image.");
00107         }
00108 
00109         if (!is_ri()) {
00110                 throw ImageFormatException("complex ri expected. Got amplitude/phase.");
00111         }
00112 
00113         int offset = is_fftodd() ? 1 : 2;
00114         EMData* dat = new EMData();
00115         int ndim = get_ndim();
00116         dat->set_size_cuda(nx-offset, ny, nz);
00117         float *d = dat->get_cuda_data();
00118         float *this_d;
00119         EMData* tmp = 0;
00120         if (preserve_input){
00121                 tmp = new EMData(*this);
00122                 this_d = tmp->get_cuda_data();
00123         } else {
00124                 this_d = get_cuda_data();
00125         }
00126         if ( ndim == 1 ) {
00127                 cuda_dd_fft_complex_to_real_nd(this_d,d, nx-offset,1,1);
00128         } else if (ndim == 2) {
00129                 cuda_dd_fft_complex_to_real_nd(this_d,d, ny,nx-offset,1);
00130         } else if (ndim == 3) {
00131                 cuda_dd_fft_complex_to_real_nd(this_d,d, nz,ny,nx-offset);
00132         } else throw ImageDimensionException("No cuda FFT support of images with dimensions exceeding 3");
00133 
00134         if (tmp != 0) delete tmp;
00135 
00136         // SCALE the inverse FFT
00137         float scale = 1.0f/static_cast<float>((dat->get_size()));
00138         dat->mult(scale); // Use of GPU should be automatic
00139 
00140         dat->set_fftpad(false);
00141         dat->set_fftodd(false);
00142         dat->set_complex(false);
00143         if(dat->get_ysize()==1 && dat->get_zsize()==1)  dat->set_complex_x(false);
00144         dat->set_ri(false);
00145         dat->gpu_update();
00146 
00147         EXITFUNC;
00148         return dat;
00149 }
00150 
00151 #endif //EMAN2_USING_CUDA
00152 
00153 EMData *EMData::do_fft() const
00154 {
00155         ENTERFUNC;
00156 
00157 #ifdef EMAN2_USING_CUDA
00158         if (gpu_operation_preferred()) {
00159                 EXITFUNC;
00160                 return do_fft_cuda();
00161         }
00162 #endif
00163 
00164         if ( is_complex() ) {
00165                 LOGERR("real image expected. Input image is complex image.");
00166                 throw ImageFormatException("real image expected. Input image is complex image.");
00167         }
00168 
00169         int nxreal = nx;
00170         int offset = 2 - nx%2;
00171         int nx2 = nx + offset;
00172         EMData* dat = copy_head();
00173         dat->set_size(nx2, ny, nz);
00174         //dat->to_zero();  // do not need it, real_to_complex will do it right anyway
00175         if (offset == 1) dat->set_fftodd(true);
00176         else             dat->set_fftodd(false);
00177 
00178         float *d = dat->get_data();
00179         //std::cout<<" do_fft "<<rdata[5]<<"  "<<d[5]<<std::endl;
00180         EMfft::real_to_complex_nd(get_data(), d, nxreal, ny, nz);
00181 
00182         dat->update();
00183         dat->set_fftpad(true);
00184         dat->set_complex(true);
00185         if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(true);
00186         dat->set_ri(true);
00187 
00188         EXITFUNC;
00189         return dat;
00190 }
00191 
00192 EMData *EMData::do_fft_inplace()
00193 {
00194         ENTERFUNC;
00195 
00196         if ( is_complex() ) {
00197                 LOGERR("real image expected. Input image is complex image.");
00198                 throw ImageFormatException("real image expected. Input image is complex image.");
00199         }
00200 
00201         size_t offset;
00202         int nxreal;
00203         get_data(); // Required call if GPU caching is being used. Otherwise harmless
00204         if (!is_fftpadded()) {
00205                 // need to extend the matrix along x
00206                 // meaning nx is the un-fftpadded size
00207                 nxreal = nx;
00208                 offset = 2 - nx%2;
00209                 if (1 == offset) set_fftodd(true);
00210                 else             set_fftodd(false);
00211                 int nxnew = nx + offset;
00212                 set_size(nxnew, ny, nz);
00213 
00214                 for (int iz = nz-1; iz >= 0; iz--) {
00215                         for (int iy = ny-1; iy >= 0; iy--) {
00216                                 for (int ix = nxreal-1; ix >= 0; ix--) {
00217                                         size_t oldxpos = ix + (iy + iz*ny)*nxreal;
00218                                         size_t newxpos = ix + (iy + iz*ny)*nxnew;
00219                                         (*this)(newxpos) = (*this)(oldxpos);
00220                                 }
00221                         }
00222                 }
00223                 set_fftpad(true);
00224         } else {
00225                 offset = is_fftodd() ? 1 : 2;
00226                 nxreal = nx - offset;
00227         }
00228         EMfft::real_to_complex_nd(rdata, rdata, nxreal, ny, nz);
00229 
00230         set_complex(true);
00231         if(ny==1 && nz==1)  set_complex_x(true);
00232         set_ri(true);
00233 
00234         update();
00235 
00236         EXITFUNC;
00237         return this;
00238 }
00239 
00240 EMData *EMData::do_ift()
00241 {
00242         ENTERFUNC;
00243 
00244 #ifdef EMAN2_USING_CUDA
00245         if (gpu_operation_preferred()) {
00246                 EXITFUNC;
00247                 return do_ift_cuda();
00248         }
00249 #endif
00250 
00251         if (!is_complex()) {
00252                 LOGERR("complex image expected. Input image is real image.");
00253                 throw ImageFormatException("complex image expected. Input image is real image.");
00254         }
00255 
00256         if (!is_ri()) {
00257                 LOGWARN("run IFT on AP data, only RI should be used. Converting.");
00258         }
00259 
00260         get_data(); // Required call if GPU caching is being used. Otherwise harmless
00261         EMData* dat = copy_head();
00262         dat->set_size(nx, ny, nz);
00263         ap2ri();
00264 
00265         float *d = dat->get_data();
00266         int ndim = get_ndim();
00267 
00268         /* Do inplace IFT on a image copy, because the complex to real transform of
00269          * nd will destroy its input array even for out-of-place transforms.
00270          */
00271         memcpy((char *) d, (char *) rdata, nx * ny * nz * sizeof(float));
00272 
00273         int offset = is_fftodd() ? 1 : 2;
00274         //cout << "Sending offset " << offset << " " << nx-offset << endl;
00275         if (ndim == 1) {
00276                 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
00277         } else {
00278                 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
00279 
00280 #ifndef CUDA_FFT
00281                 size_t row_size = (nx - offset) * sizeof(float);
00282                 for (int i = 1; i < ny * nz; i++) {
00283                         memmove((char *) &d[i * (nx - offset)], (char *) &d[i * nx], row_size);
00284                 }
00285 #endif
00286         }
00287 
00288         dat->set_size(nx - offset, ny, nz);     //remove the padding
00289 #if defined     FFTW2 || defined FFTW3  || defined CUDA_FFT //native fft and ACML already done normalization
00290         // SCALE the inverse FFT
00291         float scale = 1.0f / ((nx - offset) * ny * nz);
00292         dat->mult(scale);
00293 #endif  //FFTW2 || FFTW3
00294         dat->set_fftodd(false);
00295         dat->set_fftpad(false);
00296         dat->set_complex(false);
00297         if(dat->get_ysize()==1 && dat->get_zsize()==1)  dat->set_complex_x(false);
00298         dat->set_ri(false);
00299         dat->update();
00300 
00301 
00302         EXITFUNC;
00303         return dat;
00304 }
00305 
00306 /*
00307    FFT in place does not depad, return real x-extended image
00308    use
00309 */
00310 EMData *EMData::do_ift_inplace()
00311 {
00312         ENTERFUNC;
00313 
00314         if (!is_complex()) {
00315                 LOGERR("complex image expected. Input image is real image.");
00316                 throw ImageFormatException("complex image expected. Input image is real image.");
00317         }
00318 
00319         if (!is_ri()) {
00320                 LOGWARN("run IFT on AP data, only RI should be used. ");
00321         }
00322         ap2ri();
00323 
00324         int offset = is_fftodd() ? 1 : 2;
00325         float* data = get_data();
00326         EMfft::complex_to_real_nd(data, data, nx - offset, ny, nz);
00327 
00328 #if defined     FFTW2 || defined FFTW3 || defined CUDA_FFT      //native fft and ACML already done normalization
00329         // SCALE the inverse FFT
00330         int nxo = nx - offset;
00331         float scale = 1.0f / (nxo * ny * nz);
00332         mult(scale);
00333 #endif //FFTW2 || FFTW3
00334 
00335 #ifndef CUDA_FFT
00336         set_fftpad(true);
00337 #else
00338         set_size(nx - offset, ny, nz);
00339 #endif
00340         set_complex(false);
00341         if(ny==1 && nz==1) set_complex_x(false);
00342         set_ri(false);
00343         update();
00344 
00345 
00346         EXITFUNC;
00347         return this;
00348 }
00349 #undef rdata
00350 
00351 
00352 std::string EMData::render_ap24(int x0, int y0, int ixsize, int iysize,
00353                                                  int bpl, float scale, int mingray, int maxgray,
00354                                                  float render_min, float render_max,float gamma,int flags)
00355 {
00356         ENTERFUNC;
00357 
00358         int asrgb;
00359         int hist=(flags&2)/2;
00360         int invy=(flags&4)?1:0;
00361 
00362         if (!is_complex()) throw ImageDimensionException("complex only");
00363 
00364         if (get_ndim() != 2) {
00365                 throw ImageDimensionException("2D only");
00366         }
00367 
00368         if (is_complex()) ri2ap();
00369 
00370         if (render_max <= render_min) {
00371                 render_max = render_min + 0.01f;
00372         }
00373 
00374         if (gamma<=0) gamma=1.0;
00375 
00376         // Calculating a full floating point gamma for
00377         // each pixel in the image slows rendering unacceptably
00378         // however, applying a gamma-mapping to an 8 bit colorspace
00379         // has unaccepable accuracy. So, we oversample the 8 bit colorspace
00380         // as a 12 bit colorspace and apply the gamma mapping to that
00381         // This should produce good accuracy for gamma values
00382         // larger than 0.5 (and a high upper limit)
00383         static int smg0=0,smg1=0;       // while this destroys threadsafety in the rendering process
00384         static float sgam=0;            // it is necessary for speed when rendering large numbers of small images
00385         static unsigned char gammamap[4096];
00386         if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
00387                 for (int i=0; i<4096; i++) {
00388                         if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00389                         else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00390                 }
00391         }
00392         smg0=mingray;   // so we don't recompute the map unless something changes
00393         smg1=maxgray;
00394         sgam=gamma;
00395 
00396         if (flags&8) asrgb=4;
00397         else if (flags&1) asrgb=3;
00398         else throw ImageDimensionException("must set flag 1 or 8");
00399 
00400         std::string ret=std::string();
00401 //      ret.resize(iysize*bpl);
00402         ret.assign(iysize*bpl+hist*1024,char(mingray));
00403         unsigned char *data=(unsigned char *)ret.data();
00404         unsigned int *histd=(unsigned int *)(data+iysize*bpl);
00405         if (hist) {
00406                 for (int i=0; i<256; i++) histd[i]=0;
00407         }
00408 
00409         float rm = render_min;
00410         float inv_scale = 1.0f / scale;
00411         int ysize = iysize;
00412         int xsize = ixsize;
00413 
00414         int ymin = 0;
00415         if (iysize * inv_scale > ny) {
00416                 ymin = (int) (iysize - ny / inv_scale);
00417         }
00418 
00419         float gs = (maxgray - mingray) / (render_max - render_min);
00420         float gs2 = 4095.999f / (render_max - render_min);
00421 //      float gs2 = 1.0 / (render_max - render_min);
00422         if (render_max < render_min) {
00423                 gs = 0;
00424                 rm = FLT_MAX;
00425         }
00426 
00427         int dsx = -1;
00428         int dsy = 0;
00429         int remx = 0;
00430         int remy = 0;
00431         const int scale_n = 100000;
00432 
00433         int addi = 0;
00434         int addr = 0;
00435         if (inv_scale == floor(inv_scale)) {
00436                 dsx = (int) inv_scale;
00437                 dsy = (int) (inv_scale * nx);
00438         }
00439         else {
00440                 addi = (int) floor(inv_scale);
00441                 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00442         }
00443 
00444         int xmin = 0;
00445         if (x0 < 0) {
00446                 xmin = (int) (-x0 / inv_scale);
00447                 xsize -= (int) floor(x0 / inv_scale);
00448                 x0 = 0;
00449         }
00450 
00451         if ((xsize - xmin) * inv_scale > (nx - x0)) {
00452                 xsize = (int) ((nx - x0) / inv_scale + xmin);
00453         }
00454         int ymax = ysize - 1;
00455         if (y0 < 0) {
00456                 ymax = (int) (ysize + y0 / inv_scale - 1);
00457                 ymin += (int) floor(y0 / inv_scale);
00458                 y0 = 0;
00459         }
00460 
00461         if (xmin < 0) xmin = 0;
00462         if (ymin < 0) ymin = 0;
00463         if (xsize > ixsize) xsize = ixsize;
00464         if (ymax > iysize) ymax = iysize;
00465 
00466         int lmax = nx * ny - 1;
00467 
00468         int mid=nx*ny/2;
00469         float* image_data = get_data();
00470         if (dsx != -1) {
00471                 int l = y0 * nx;
00472                 for (int j = ymax; j >= ymin; j--) {
00473                         int ll = x0;
00474                         for (int i = xmin; i < xsize; i++) {
00475                                 if (l + ll > lmax || ll >= nx - 2) break;
00476 
00477                                 int k = 0;
00478                                 unsigned char p;
00479                                 int ph;
00480                                 if (ll >= nx / 2) {
00481                                         if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
00482                                         else k = 2 * (ll - nx / 2) + l + 2 + nx;
00483                                         ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
00484                                 }
00485                                 else {
00486                                         k = nx * ny - (l + 2 * ll) - 2;
00487                                         ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;        // complex phase as integer 0-767
00488                                 }
00489                                 if (k>=mid) k-=mid;             // These 2 lines handle the Fourier origin being in the corner, not the middle
00490                                 else k+=mid;
00491                                 float t = image_data[k];
00492                                 if (t <= rm)  p = mingray;
00493                                 else if (t >= render_max) p = maxgray;
00494                                 else if (gamma!=1.0) {
00495                                         k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00496                                         p = gammamap[k];                                        // apply gamma using precomputed gamma map
00497                                 }
00498                                 else {
00499                                         p = (unsigned char) (gs * (t - render_min));
00500                                         p += mingray;
00501                                 }
00502                                 if (ph<256) {
00503                                         data[i * asrgb + j * bpl] = p*(255-ph)/256;
00504                                         data[i * asrgb + j * bpl+1] = p*ph/256;
00505                                         data[i * asrgb + j * bpl+2] = 0;
00506                                 }
00507                                 else if (ph<512) {
00508                                         data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00509                                         data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00510                                         data[i * asrgb + j * bpl] = 0;
00511                                 }
00512                                 else {
00513                                         data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00514                                         data[i * asrgb + j * bpl] = p*(ph-512)/256;
00515                                         data[i * asrgb + j * bpl+1] = 0;
00516                                 }
00517                                 if (hist) histd[p]++;
00518                                 ll += dsx;
00519                         }
00520                         l += dsy;
00521                 }
00522         }
00523         else {
00524                 remy = 10;
00525                 int l = y0 * nx;
00526                 for (int j = ymax; j >= ymin; j--) {
00527                         int br = l;
00528                         remx = 10;
00529                         int ll = x0;
00530                         for (int i = xmin; i < xsize - 1; i++) {
00531                                 if (l + ll > lmax || ll >= nx - 2) {
00532                                         break;
00533                                 }
00534                                 int k = 0;
00535                                 unsigned char p;
00536                                 int ph;
00537                                 if (ll >= nx / 2) {
00538                                         if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
00539                                         else k = 2 * (ll - nx / 2) + l + 2 + nx;
00540                                         ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
00541                                 }
00542                                 else {
00543                                         k = nx * ny - (l + 2 * ll) - 2;
00544                                         ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;        // complex phase as integer 0-767
00545                                 }
00546                                 if (k>=mid) k-=mid;             // These 2 lines handle the Fourier origin being in the corner, not the middle
00547                                 else k+=mid;
00548 
00549                                 float t = image_data[k];
00550                                 if (t <= rm)
00551                                         p = mingray;
00552                                 else if (t >= render_max) {
00553                                         p = maxgray;
00554                                 }
00555                                 else if (gamma!=1.0) {
00556                                         k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00557                                         p = gammamap[k];                                        // apply gamma using precomputed gamma map
00558                                 }
00559                                 else {
00560                                         p = (unsigned char) (gs * (t - render_min));
00561                                         p += mingray;
00562                                 }
00563                                 if (ph<256) {
00564                                         data[i * asrgb + j * bpl] = p*(255-ph)/256;
00565                                         data[i * asrgb + j * bpl+1] = p*ph/256;
00566                                         data[i * asrgb + j * bpl+2] = 0;
00567                                 }
00568                                 else if (ph<512) {
00569                                         data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00570                                         data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00571                                         data[i * asrgb + j * bpl] = 0;
00572                                 }
00573                                 else {
00574                                         data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00575                                         data[i * asrgb + j * bpl] = p*(ph-512)/256;
00576                                         data[i * asrgb + j * bpl+1] = 0;
00577                                 }
00578                                 if (hist) histd[p]++;
00579                                 ll += addi;
00580                                 remx += addr;
00581                                 if (remx > scale_n) {
00582                                         remx -= scale_n;
00583                                         ll++;
00584                                 }
00585                         }
00586                         l = br + addi * nx;
00587                         remy += addr;
00588                         if (remy > scale_n) {
00589                                 remy -= scale_n;
00590                                 l += nx;
00591                         }
00592                 }
00593         }
00594 
00595         // this replicates r -> g,b
00596         if (asrgb==4) {
00597                 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00598                         for (int i=xmin; i<xsize*4; i+=4) {
00599                                 data[i+j+3]=255;
00600                         }
00601                 }
00602         }
00603 
00604         EXITFUNC;
00605 
00606         // ok, ok, not the most efficient place to do this, but it works
00607         if (invy) {
00608                 int x,y;
00609                 char swp;
00610                 for (y=0; y<iysize/2; y++) {
00611                         for (x=0; x<ixsize; x++) {
00612                                 swp=ret[y*bpl+x];
00613                                 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
00614                                 ret[(iysize-y-1)*bpl+x]=swp;
00615                         }
00616                 }
00617         }
00618 
00619     //  return PyString_FromStringAndSize((const char*) data,iysize*bpl);
00620         return ret;
00621 }
00622 
00623 
00624 void EMData::render_amp24( int x0, int y0, int ixsize, int iysize,
00625                                                   int bpl, float scale, int mingray, int maxgray,
00626                                                   float render_min, float render_max, void *ref,
00627                                                   void cmap(void *, int coord, unsigned char *tri))
00628 {
00629         ENTERFUNC;
00630 
00631         if (get_ndim() != 2) {
00632                 throw ImageDimensionException("2D only");
00633         }
00634 
00635         if (is_complex()) {
00636                 ri2ap();
00637         }
00638 
00639         if (render_max <= render_min) {
00640                 render_max = render_min + 0.01f;
00641         }
00642 
00643         std::string ret=std::string();
00644         ret.resize(iysize*bpl);
00645         unsigned char *data=(unsigned char *)ret.data();
00646 
00647         float rm = render_min;
00648         float inv_scale = 1.0f / scale;
00649         int ysize = iysize;
00650         int xsize = ixsize;
00651         const int scale_n = 100000;
00652 
00653         int ymin = 0;
00654         if ( iysize * inv_scale > ny) {
00655                 ymin = (int) (iysize - ny / inv_scale);
00656         }
00657         float gs = (maxgray - mingray) / (render_max - render_min);
00658         if (render_max < render_min) {
00659                 gs = 0;
00660                 rm = FLT_MAX;
00661         }
00662         int dsx = -1;
00663         int dsy = 0;
00664         if (inv_scale == floor(inv_scale)) {
00665                 dsx = (int) inv_scale;
00666                 dsy = (int) (inv_scale * nx);
00667         }
00668         int addi = 0;
00669         int addr = 0;
00670 
00671         if (dsx == -1) {
00672                 addi = (int) floor(inv_scale);
00673                 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00674         }
00675 
00676         int remx = 0;
00677         int remy = 0;
00678         int xmin = 0;
00679         if (x0 < 0) {
00680                 xmin = (int) (-x0 / inv_scale);
00681                 xsize -= (int) floor(x0 / inv_scale);
00682                 x0 = 0;
00683         }
00684 
00685         if ((xsize - xmin) * inv_scale > (nx - x0)) {
00686                 xsize = (int) ((nx - x0) / inv_scale + xmin);
00687         }
00688         int ymax = ysize - 1;
00689         if (y0 < 0) {
00690                 ymax = (int) (ysize + y0 / inv_scale - 1);
00691                 ymin += (int) floor(y0 / inv_scale);
00692                 y0 = 0;
00693         }
00694 
00695 
00696         if (xmin < 0) {
00697                 xmin = 0;
00698         }
00699 
00700         if (ymin < 0) {
00701                 ymin = 0;
00702         }
00703         if (xsize > ixsize) {
00704                 xsize = ixsize;
00705         }
00706         if (ymax > iysize) {
00707                 ymax = iysize;
00708         }
00709 
00710         int lmax = nx * ny - 1;
00711         unsigned char tri[3];
00712         float* image_data = get_data();
00713         if (is_complex()) {
00714                 if (dsx != -1) {
00715                         int l = y0 * nx;
00716                         for (int j = ymax; j >= ymin; j--) {
00717                                 int ll = x0;
00718                                 for (int i = xmin; i < xsize; i++, ll += dsx) {
00719                                         if (l + ll > lmax || ll >= nx - 2) {
00720                                                 break;
00721                                         }
00722                                         int kk = 0;
00723                                         if (ll >= nx / 2) {
00724                                                 if (l >= (ny - inv_scale) * nx) {
00725                                                         kk = 2 * (ll - nx / 2) + 2;
00726                                                 }
00727                                                 else {
00728                                                         kk = 2 * (ll - nx / 2) + l + 2 + nx;
00729                                                 }
00730                                         }
00731                                         else {
00732                                                 kk = nx * ny - (l + 2 * ll) - 2;
00733                                         }
00734                                         int k = 0;
00735                                         float t = image_data[kk];
00736                                         if (t <= rm) {
00737                                                 k = mingray;
00738                                         }
00739                                         else if (t >= render_max) {
00740                                                 k = maxgray;
00741                                         }
00742                                         else {
00743                                                 k = (int) (gs * (t - render_min));
00744                                                 k += mingray;
00745                                         }
00746                                         tri[0] = static_cast < unsigned char >(k);
00747                                         cmap(ref, kk, tri);
00748                                         data[i * 3 + j * bpl] = tri[0];
00749                                         data[i * 3 + 1 + j * bpl] = tri[1];
00750                                         data[i * 3 + 2 + j * bpl] = tri[2];
00751                                 }
00752                                 l += dsy;
00753                         }
00754                 }
00755                 else {
00756                         remy = 10;
00757                         for (int j = ymax, l = y0 * nx; j >= ymin; j--) {
00758                                 int br = l;
00759                                 remx = 10;
00760                                 for (int i = xmin, ll = x0; i < xsize - 1; i++) {
00761                                         if (l + ll > lmax || ll >= nx - 2) {
00762                                                 break;
00763                                         }
00764                                         int kk = 0;
00765                                         if (ll >= nx / 2) {
00766                                                 if (l >= (ny * nx - nx)) {
00767                                                         kk = 2 * (ll - nx / 2) + 2;
00768                                                 }
00769                                                 else {
00770                                                         kk = 2 * (ll - nx / 2) + l + 2 + nx;
00771                                                 }
00772                                         }
00773                                         else {
00774                                                 kk = nx * ny - (l + 2 * ll) - 2;
00775                                         }
00776                                         int k = 0;
00777                                         float t = image_data[kk];
00778                                         if (t <= rm) {
00779                                                 k = mingray;
00780                                         }
00781                                         else if (t >= render_max) {
00782                                                 k = maxgray;
00783                                         }
00784                                         else {
00785                                                 k = (int) (gs * (t - render_min));
00786                                                 k += mingray;
00787                                         }
00788                                         tri[0] = static_cast < unsigned char >(k);
00789                                         cmap(ref, kk, tri);
00790                                         data[i * 3 + j * bpl] = tri[0];
00791                                         data[i * 3 + 1 + j * bpl] = tri[1];
00792                                         data[i * 3 + 2 + j * bpl] = tri[2];
00793                                         ll += addi;
00794                                         remx += addr;
00795                                         if (remx > scale_n) {
00796                                                 remx -= scale_n;
00797                                                 ll++;
00798                                         }
00799                                 }
00800                                 l = br + addi * nx;
00801                                 remy += addr;
00802                                 if (remy > scale_n) {
00803                                         remy -= scale_n;
00804                                         l += nx;
00805                                 }
00806                         }
00807                 }
00808         }
00809         else {
00810                 if (dsx != -1) {
00811                         for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
00812                                 int br = l;
00813                                 for (int i = xmin; i < xsize; i++, l += dsx) {
00814                                         if (l > lmax) {
00815                                                 break;
00816                                         }
00817                                         float t = image_data[l];
00818                                         int k = 0;
00819                                         if (t <= rm) {
00820                                                 k = mingray;
00821                                         }
00822                                         else if (t >= render_max) {
00823                                                 k = maxgray;
00824                                         }
00825                                         else {
00826                                                 k = (int) (gs * (t - render_min));
00827                                                 k += mingray;
00828                                         }
00829                                         tri[0] = static_cast < unsigned char >(k);
00830                                         cmap(ref, l, tri);
00831                                         data[i * 3 + j * bpl] = tri[0];
00832                                         data[i * 3 + 1 + j * bpl] = tri[1];
00833                                         data[i * 3 + 2 + j * bpl] = tri[2];
00834                                 }
00835                                 l = br + dsy;
00836                         }
00837                 }
00838                 else {
00839                         remy = 10;
00840                         for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
00841                                 int br = l;
00842                                 remx = 10;
00843                                 for (int i = xmin; i < xsize; i++) {
00844                                         if (l > lmax) {
00845                                                 break;
00846                                         }
00847                                         float t = image_data[l];
00848                                         int k = 0;
00849                                         if (t <= rm) {
00850                                                 k = mingray;
00851                                         }
00852                                         else if (t >= render_max) {
00853                                                 k = maxgray;
00854                                         }
00855                                         else {
00856                                                 k = (int) (gs * (t - render_min));
00857                                                 k += mingray;
00858                                         }
00859                                         tri[0] = static_cast < unsigned char >(k);
00860                                         cmap(ref, l, tri);
00861                                         data[i * 3 + j * bpl] = tri[0];
00862                                         data[i * 3 + 1 + j * bpl] = tri[1];
00863                                         data[i * 3 + 2 + j * bpl] = tri[2];
00864                                         l += addi;
00865                                         remx += addr;
00866                                         if (remx > scale_n) {
00867                                                 remx -= scale_n;
00868                                                 l++;
00869                                         }
00870                                 }
00871                                 l = br + addi * nx;
00872                                 remy += addr;
00873                                 if (remy > scale_n) {
00874                                         remy -= scale_n;
00875                                         l += nx;
00876                                 }
00877                         }
00878                 }
00879         }
00880 
00881         EXITFUNC;
00882 }
00883 
00884 void EMData::ap2ri()
00885 {
00886         ENTERFUNC;
00887 
00888         if (!is_complex() || is_ri()) {
00889                 return;
00890         }
00891 
00892 #ifdef EMAN2_USING_CUDA
00893         if (gpu_operation_preferred()) {
00894                 EMDataForCuda tmp = get_data_struct_for_cuda();
00895                 emdata_ap2ri(&tmp);
00896                 set_ri(true);
00897                 gpu_update();
00898                 EXITFUNC;
00899                 return;
00900         }
00901 #endif
00902 
00903         Util::ap2ri(get_data(), nx * ny * nz);
00904         set_ri(true);
00905         update();
00906         EXITFUNC;
00907 }
00908 
00909 void EMData::ri2inten()
00910 {
00911         ENTERFUNC;
00912 
00913         if (!is_complex()) return;
00914         if (!is_ri()) ap2ri();
00915 
00916 #ifdef EMAN2_USING_CUDA
00917         if (gpu_operation_preferred()) {
00918                 EMDataForCuda tmp = get_data_struct_for_cuda();
00919                 emdata_ri2inten(&tmp);
00920                 set_attr("is_intensity", int(1));
00921                 gpu_update();
00922                 EXITFUNC;
00923                 return;
00924         }
00925 #endif
00926 
00927         float * data = get_data();
00928         size_t size = nx * ny * nz;
00929         for (size_t i = 0; i < size; i += 2) {
00930                 data[i]=data[i]*data[i]+data[i+1]*data[i+1];
00931                 data[i+1]=0;
00932         }
00933 
00934         set_attr("is_intensity", int(1));
00935         update();
00936         EXITFUNC;
00937 }
00938 
00939 
00940 void EMData::ri2ap()
00941 {
00942         ENTERFUNC;
00943 
00944         if (!is_complex() || !is_ri()) {
00945                 return;
00946         }
00947 #ifdef EMAN2_USING_CUDA
00948         if (gpu_operation_preferred()) {
00949                 EMDataForCuda tmp = get_data_struct_for_cuda();
00950                 emdata_ri2ap(&tmp);
00951                 set_ri(false);
00952                 gpu_update();
00953                 EXITFUNC;
00954                 return;
00955         }
00956 #endif
00957 
00958         float * data = get_data();
00959 
00960         size_t size = nx * ny * nz;
00961         for (size_t i = 0; i < size; i += 2) {
00962 #ifdef  _WIN32
00963                 float f = (float)_hypot(data[i], data[i + 1]);
00964 #else
00965                 float f = (float)hypot(data[i], data[i + 1]);
00966 #endif
00967                 if (data[i] == 0 && data[i + 1] == 0) {
00968                         data[i + 1] = 0;
00969                 }
00970                 else {
00971                         data[i + 1] = atan2(data[i + 1], data[i]);
00972                 }
00973                 data[i] = f;
00974         }
00975 
00976         set_ri(false);
00977         update();
00978         EXITFUNC;
00979 }
00980 
00981 
00982 float calc_bessel(const int n, const float& x) {
00983         gsl_sf_result result;
00984 //      int success =
00985         gsl_sf_bessel_Jn_e(n,(double)x, &result);
00986         return (float)result.val;
00987 }
00988 
00989 EMData*   EMData::bispecRotTransInvN(int N, int NK)
00990 {
00991 
00992         int EndP = this -> get_xsize(); // length(fTrueVec);
00993         int Mid  = (int) ((1+EndP)/2);
00994         int End = 2*Mid-1;
00995 
00996         int CountxyMax = End*End;
00997 
00998         int   *SortfkInds       = new    int[CountxyMax];
00999         int   *kVecX            = new    int[CountxyMax];
01000         int   *kVecY            = new    int[CountxyMax];
01001         float *fkVecR           = new  float[CountxyMax];
01002         float *fkVecI           = new  float[CountxyMax];
01003         float *absD1fkVec       = new  float[CountxyMax];
01004         float *absD1fkVecSorted = new  float[CountxyMax];
01005 
01006         float *jxjyatan2         = new  float[End*End]; //  jxjyatan2[jy*End + jx]  = atan2(jy+1-Mid , jx +1 -Mid);
01007 
01008         EMData * ThisCopy = new EMData(End,End);
01009 
01010         for (int jx=0; jx <End ; jx++) {
01011                 for (int jy=0; jy <End ; jy++) {
01012                         float ValNow = this -> get_value_at(jx,jy);
01013                         ThisCopy -> set_value_at(jx,jy,ValNow);
01014 //              cout<< " jxM= " << jx+1<<" jyM= " << jy+1<< "ValNow" << ValNow << endl; //    Works
01015         }}
01016 
01017 
01018         EMData* fk = ThisCopy -> do_fft();
01019         fk          ->process_inplace("xform.fourierorigin.tocenter");
01020 
01021 //      EMData* fk
01022         EMData* fkRCopy = new EMData(End,End);
01023         EMData* fkICopy = new EMData(End,End);
01024         EMData* fkCopy  = new EMData(End,End);
01025 
01026 
01027         for  (int jCount= 0; jCount<End*End; jCount++) {
01028 //              jCount = jy*End + jx;
01029                 int jx             = jCount%End ;
01030                 int jy             = (jCount-jx)/End ;
01031                 jxjyatan2[jCount]  = atan2((float)(jy+1-Mid) , (float)(jx +1-Mid));
01032         }
01033 
01034 
01035         for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) { // kEx twice the value of the Fourier
01036                                                 // x variable: EMAN index for real, imag
01037                 int kx    = kEx/2;              // kx  is  the value of the Fourier variable
01038                 int kIx   = kx+Mid-1; // This is the value of the index for a matlab image (-1)
01039                 int kCx   =  -kx ;
01040                 int kCIx  = kCx+ Mid-1 ;
01041                 for (int kEy= 0 ; kEy<End; kEy++) { // This is the value of the EMAN index
01042                         int kIy              =  kEy       ; //  This is the value of the index for a matlab image (-1)
01043                         int ky               =  kEy+1-Mid; // (kEy+ Mid-1)%End - Mid+1 ;  // This is the actual value of the Fourier variable
01044                         float realVal        =  fk -> get_value_at(kEx  ,kEy) ;
01045                         float imagVal        =  fk -> get_value_at(kEx+1,kEy) ;
01046                         float absVal         =  ::sqrt(realVal*realVal+imagVal*imagVal);
01047                         float fkAng          =  atan2(imagVal,realVal);
01048 
01049                         float NewRealVal   ;
01050                         float NewImagVal   ;
01051                         float AngMatlab    ;
01052 
01053                         if (kIx==Mid-1) {
01054 //                              AngMatlab = -fkAng - 2.*M_PI*(kIy+ 1-Mid)*(Mid)/End;
01055 //                      cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]="  << fkVecI[i] <<"  angle[i]= "  << AngMatlab << endl;
01056                         }
01057 
01058                         if (kIx>Mid-1){
01059 //                      cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]="  << fkVecI[i] <<"  angle[i]= "  << AngMatlab << endl;
01060                         }
01061 
01062                         AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
01063                         NewRealVal  =   absVal*cos(AngMatlab);
01064                         NewImagVal  =   absVal*sin(AngMatlab);
01065 
01066 
01067                         fkVecR[kIy+kIx *End] =  NewRealVal ;
01068                         fkVecR[kIy+kCIx*End] =  NewRealVal ;
01069                         fkVecI[kIy+kIx *End] =  NewImagVal ;
01070                         fkVecI[kIy+kCIx*End] = -NewImagVal ;
01071                         absD1fkVec[kIy + kIx  *End] = absVal;
01072                         absD1fkVec[kIy + kCIx *End] = absVal;
01073                         kVecX[kIy+kIx  *End] =  kx      ;
01074                         kVecX[kIy+kCIx *End] =  kCx    ;
01075                         kVecY[kIy+kIx  *End] =  ky     ;
01076                         kVecY[kIy+kCIx *End] =  ky     ;
01077 //                      printf("kx=%d,ky=%d,tempVal =%f+ i %4.2f \n",kx,ky,realVal,imagVal );
01078 //                      cout << "kx = " << kx << "; ky = "<< ky << "; val is" << realVal<<"+ i "<<imagVal<< endl;
01079 
01080 //                      cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; fkAng*9/ 2pi is " << fkAng*9/2/M_PI<<  endl;
01081 //                      cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; absval is " << absVal<<  "; realval is " << NewRealVal<< "; imagval is " << NewImagVal<< endl;
01082                         fkCopy  -> set_value_at(kIx ,kIy, absVal);
01083                         fkCopy  -> set_value_at(kCIx,kIy, absVal);
01084                         fkRCopy -> set_value_at(kIx, kIy, NewRealVal);
01085                         fkRCopy -> set_value_at(kCIx,kIy, NewRealVal);
01086                         fkICopy -> set_value_at(kIx, kIy, NewImagVal);
01087                         fkICopy -> set_value_at(kCIx,kIy,-NewImagVal);
01088 
01089                 }
01090         }
01091         system("rm -f fkCopy.???");
01092         system("rm -f fk?Copy.???");
01093         fkCopy  -> write_image("fkCopy.img");
01094         fkRCopy -> write_image("fkRCopy.img");
01095         fkICopy -> write_image("fkICopy.img");
01096 
01097         cout << "Starting the sort "<< endl;
01098 
01099         vector< pair<float, int> > absInds;
01100         for(int i  = 0; i < CountxyMax; ++i ) {
01101                 pair<float,int> p;
01102                 p = make_pair(absD1fkVec[i],i); // p = make_pair(rand(),i);
01103                 absInds.push_back( p);
01104         }
01105 
01106         std::sort(absInds.begin(),absInds.end());
01107 
01108         for(int i  = 0; i < CountxyMax; ++i ) {
01109                 pair<float,int> p   ;
01110                 p = absInds[i]         ;
01111                 absD1fkVecSorted[CountxyMax-1-i] =  p.first ;
01112                 SortfkInds[CountxyMax-1-i]       =  p.second ;
01113         }
01114 
01115         cout << "Ending the sort "<< endl;
01116 
01117 // float AngsMat[] ={2.8448, -0.3677,-0.2801,-1.0494,-1.7836,-2.5179, 2.9959, 3.0835,-0.1290,-0.8876,2.1829, 2.2705,1.5011,0.7669,0.0327,-0.7366,-0.6489,2.4215,-1.6029,1.4676,1.5552,0.7859,0.0517,-0.6825,-1.4518,-1.3642,1.7063,-1.7845,1.2859,1.3736,0.6043,-0.1299,-0.8642,-1.6335,-1.5459,1.5247,-1.6546,1.4159,1.5036,0.7342,0,-0.7342,-1.5036,-1.4159,1.6546,-1.5247,1.5459,1.6335,0.8642,0.1299,-0.6043,-1.3736,-1.286,1.7846,-1.7063,1.3642,1.4519,0.6825,-0.0517,-0.7859,-1.5553,-1.4676,1.6029,-2.4216,0.649,0.7366,-0.0327,-0.767,-1.5012,-2.2705,-2.1829,0.8877,0.1291,-3.0836,-2.9959,2.5179,1.7837,1.0495,0.2801,0.3677,-2.8449};
01118 
01119 
01120         for(int i  = 0; i < CountxyMax; ++i ) {  // creates a new fkVec
01121                 int Si  = SortfkInds[i];
01122                 int kIx = (int)  Si/End;  kIx = (int)  i/End; // i = kIx*End+kIy
01123                 int kIy = Si  - kIx*End;  kIy = i  - kIx*End;
01124                 int iC = (End-1-kIx)*End + (End-1-kIy);
01125 //              if (i<30) { cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " valAft=" << absD1fkVecSorted[i]<< " valBef="  <<     absD1fkVec[Si] << "  SortfkInds = " << Si <<endl; }// This worked
01126 //              cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]="  << fkVecI[i] <<"  angle[i]= "  << fkAng << endl;
01127         }
01128         cout<< "Ratio of Last Amplitude to First Amplitude= " << absD1fkVecSorted[NK] /absD1fkVecSorted[0]  << endl;
01129 
01130 //      pause;
01131 
01132         for(int i  = 0; i < NK; ++i ) { // Prints out the new fkVec ,  CountxyMax
01133                 int Si= SortfkInds[i];
01134                 int kIx = (int)  Si/End; // i = kIx*End+kIy
01135                 int kIy = Si  - kIx*End;
01136  //             cout << " kIxM= " << kIx+1 << " kIyM=" << kIy+1 << " fkVecAbs=" << ::sqrt(fkVecR[Si]*fkVecR[Si] +  fkVecI[Si]* fkVecI[Si]) << " fkVecAbs=" << absD1fkVecSorted[i] << " kx= " << kVecX[Si] <<  " ky=" << kVecY[Si] <<  endl;
01137         }
01138 
01139 //       angEMAN+angMat+angDiff    =0  mod 2 pi
01140 
01141 //      angDiff=  2*pi*((-4):4)*(Mid)/End; angEMAN+angMat+angDiff= integer*2 *pi
01142 //              [  absD1fkVecSorted, SortfkInds] =sort( absD1fkVec,'descend') ;
01143 //      Util::sort_mat(&absD1fkVec[0],&absD1fkVec[Countxy],&SortfkInds[0],&SortfkInds[Countxy]);
01144 
01145 
01146 //      Let radial sampling be 0:0.5:(Mid-1)
01147 
01148  //     int NK=  min(12,CountxyMax) ;
01149 
01150 
01151 
01152         cout << "NK = " << NK << endl;
01153         float frR= 3.0/4.0;
01154         int LradRange= (int) (floor(Mid/frR)) ;
01155 
01156         float *radRange = new float[LradRange]; //= 0:.75:(Mid-1);
01157         radRange[0]=0;
01158         for (int irad=1; irad < LradRange; irad++){
01159                         radRange[irad] = radRange[irad-1] + frR; }
01160 
01161 
01162 
01163          // should equal to (2*Mid-1)
01164         cout << "Starting the calculation of invariants for N= " << N << endl;
01165 
01166 /*      int NMax=5;            */
01167 
01168         EMData* RotTransInv = new EMData();
01169         RotTransInv -> set_size(LradRange,LradRange);
01170 
01171 
01172 //      float  *RotTransInv       = new float[LradRange*LradRange ] ;
01173 //      float  *RotTransInvN      = new float[LradRange*LradRange*(NMax+1) ] ;
01174 
01175 //      for (int N=0 ; N<NMax; N++) {
01176 
01177         for (int jr1=0; jr1 < LradRange ; jr1++ ) { // LradRange
01178                 float r1= radRange[jr1];
01179 //              cout << "Pre jr2 "<< endl;
01180                 for (int jr2=0;  jr2<LradRange;  jr2++ ) { //LradRange
01181                         float r2= radRange[jr2];
01182                         float RotTransInvTemp=0;
01183                         for (int jCountkxy =0; jCountkxy<NK; jCountkxy++){
01184                                 int Countkxy =SortfkInds[jCountkxy] ;   // 1: CountxyMax
01185                                 int kx = kVecX[Countkxy] ;
01186                                 int ky = kVecY[Countkxy] ;
01187                                 float k2 = (float)(kx*kx+ky*ky);
01188                                 if (k2==0) { continue;}
01189                                 float phiK =0;
01190                                 if (k2>0) phiK= jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1];  phiK= atan2((float)ky,(float)kx);
01191 
01192                                 float fkR     = fkVecR[Countkxy] ;
01193                                 float fkI     = fkVecI[Countkxy]  ;
01194 /*                              printf("jCountkxy=%d, Countkxy=%d,absD1fkVec(Countkxy)=%f,\t\t kx=%d, ky=%d \n", jCountkxy, Countkxy, absD1fkVec[Countkxy], kx, ky);*/
01195 
01196                                 for (int jCountqxy =0; jCountqxy<NK; jCountqxy++){
01197                                         int Countqxy =SortfkInds[jCountqxy] ;   // Countqxy is the index for absD1fkVec
01198                                         int qx   = kVecX[Countqxy] ;
01199                                         int qy   = kVecY[Countqxy] ;
01200                                         int q2   = qx*qx+qy*qy;
01201                                         if (q2==0) {continue;}
01202                                         float phiQ =0;
01203                                         if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1];   phiQ=atan2((float)qy,(float)qx);
01204                                         float fqR     = fkVecR[Countqxy]  ;
01205                                         float fqI     = fkVecI[Countqxy]  ;
01206                                         int kCx  = (-kx-qx);
01207                                         int kCy  = (-ky-qy);
01208                                         int kCIx = ((kCx+Mid+2*End)%End);// labels of the image in C
01209                                         int kCIy = ((kCy+Mid+2*End)%End);
01210                                         kCx  = kCIx-Mid; // correct
01211                                         kCy  = kCIy-Mid; // correct
01212                                         int CountCxy = kCIx*End+kCIy;
01213                                         float fCR     = fkVecR[CountCxy];
01214                                         float fCI     = fkVecI[CountCxy];
01215                                         if (jr1+jr2==-1) {
01216                                         printf("jCountqxy=%d , Countqxy=%d, absD1fkVec(Countqxy)=%f,qx=%d, qy=%d \n", jCountqxy, Countqxy, absD1fkVec[Countqxy],qx, qy);
01217                                         printf(" CountCxy=%d,absD1fkVec[CountCxy]=%f,  kCx=%d,     kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
01218                                         }
01219                                         for (int p=0; p<NK; p++){
01220 //                                              printf("p=%d, SortfkInds[p]=%d, CountCxy =%d \n", p,SortfkInds[p], CountCxy);
01221                                                 if (SortfkInds[p]==CountCxy){
01222                                                         float Arg1 = 2.0f*M_PI*r1*::sqrt((float) q2)/End;
01223                                                         float Arg2 = 2.0f*M_PI*r2*::sqrt((float) k2)/End;
01224 //                                                      printf("Arg1=%4.2f, Arg2=%4.2f,  \n",Arg1, Arg2 );
01225 //                                                      if (Arg1+ Arg2<15) {
01226                                                                 float bispectemp  = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR  +fqR*fCI))
01227                                                                 * cos(N*(phiK-phiQ+M_PI));
01228                                                                 bispectemp  -= (fkR*(fqR*fCI + fqI*fCR) +fkI*(fqR*fCR - fqI*fCI))
01229                                                                 * sin(N*(phiK-phiQ+M_PI));
01230                                                                 float bess1 = calc_bessel(N, Arg1 );
01231                                                                 float bess2 = calc_bessel(N, Arg2 );
01232 //                      printf("fkr=%4.2f, fqr=%4.2f, bess1=%4.2f,bess2=%4.2f \n",fkR, fqR, bess1, bess2);
01233 /*                      printf("p =%d, SortfkInds[p]=%d, CountCxy=%d, Arg1 =%4.2f, bess1=%4.2f,  \n",
01234                                 p, SortfkInds[p],CountCxy, Arg1, bess1);*/
01235                                                                 RotTransInvTemp   = RotTransInvTemp  + bispectemp  * bess1*bess2 ;
01236 //                                                      }
01237                                                 }
01238                                         }
01239                                 } // jCountqxy
01240                         } // jCountkxy
01241                         RotTransInv -> set_value_at(jr1,jr2, RotTransInvTemp)   ;
01242 /*              RotTransInvN[jr1 + LradRange*jr2+LradRange*LradRange*N] = RotTransInvTemp  ;*/
01243                 } //jr2
01244         } //jr1
01245 // }//N
01246 
01247         return  RotTransInv ;
01248 
01249 
01250 }
01251 
01252 
01253 
01254 /*
01255 // find example
01256 #include <iostream>
01257 #include <algorithm>
01258 #include <vector>
01259 using namespace std;
01260 
01261 int main () {
01262   int myints[] = { 10, 20, 30 ,40 };
01263   int * p;
01264 
01265   // pointer to array element:
01266   p = find(myints,myints+4,30);
01267   ++p;
01268   cout << "The element following 30 is " << *p << endl;
01269 
01270   vector<int> myvector (myints,myints+4);
01271   vector<int>::iterator it;
01272 
01273   // iterator to vector element:
01274   it = find (myvector.begin(), myvector.end(), 30);
01275   ++it;
01276   cout << "The element following 30 is " << *it << endl;
01277 
01278   return 0;
01279 }*/
01280 
01281 EMData*   EMData::bispecRotTransInvDirect(int type)
01282 {
01283 
01284         int EndP = this -> get_xsize(); // length(fTrueVec);
01285         int Mid  = (int) ((1+EndP)/2);
01286         int End = 2*Mid-1;
01287 
01288         int CountxyMax = End*End;
01289 
01290         int   *SortfkInds       = new    int[CountxyMax];
01291         int   *kVecX            = new    int[CountxyMax];
01292         int   *kVecY            = new    int[CountxyMax];
01293         float *fkVecR           = new  float[CountxyMax];
01294         float *fkVecI           = new  float[CountxyMax];
01295         float *absD1fkVec       = new  float[CountxyMax];
01296         float *absD1fkVecSorted = new  float[CountxyMax];
01297 
01298 
01299         float *jxjyatan2         = new  float[End*End];
01300 
01301 
01302         EMData * ThisCopy = new EMData(End,End);
01303 
01304         for (int jx=0; jx <End ; jx++) {  // create jxjyatan2
01305                 for (int jy=0; jy <End ; jy++) {
01306                         float ValNow = this -> get_value_at(jx,jy);
01307                         ThisCopy -> set_value_at(jx,jy,ValNow);
01308                         jxjyatan2[jy*End + jx]  = atan2((float)(jy+1-Mid) , (float)(jx +1 -Mid));
01309 //              cout<< " jxM= " << jx+1<<" jyM= " << jy+1<< "ValNow" << ValNow << endl; //    Works
01310         }}
01311 
01312 
01313         EMData* fk = ThisCopy -> do_fft();
01314         fk          ->process_inplace("xform.fourierorigin.tocenter");
01315 
01316 //      Create kVecX , kVecy etc
01317 
01318         for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) { // kEx twice the value of the Fourier
01319                                                 // x variable: EMAN index for real, imag
01320                 int kx    = kEx/2;              // kx  is  the value of the Fourier variable
01321                 int kIx   = kx+Mid-1; // This is the value of the index for a matlab image (-1)
01322                 int kCx   = -kx ;
01323                 int kCIx  = kCx+ Mid-1 ;
01324                 for (int kEy= 0 ; kEy<End; kEy++) { // This is the value of the EMAN index
01325                         int kIy              =  kEy       ; //  This is the value of the index for a matlab image (-1)
01326                         int ky               =  kEy+1-Mid; // (kEy+ Mid-1)%End - Mid+1 ;  // This is the actual value of the Fourier variable
01327                         float realVal        =  fk -> get_value_at(kEx  ,kEy) ;
01328                         float imagVal        =  fk -> get_value_at(kEx+1,kEy) ;
01329                         float absVal         =  ::sqrt(realVal*realVal+imagVal*imagVal);
01330                         float fkAng          =  atan2(imagVal,realVal);
01331 
01332                         float NewRealVal   ;
01333                         float NewImagVal   ;
01334                         float AngMatlab    ;
01335 
01336                         if (kIx==Mid-1) {
01337 //                              AngMatlab = -fkAng - 2.*M_PI*(kIy+ 1-Mid)*(Mid)/End;
01338                         }
01339 
01340                         if (kIx>Mid-1){
01341 //                      cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]="  << fkVecI[i] <<"  angle[i]= "  << AngMatlab << endl;
01342                         }
01343 
01344                         AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
01345                         NewRealVal  =   absVal*cos(AngMatlab);
01346                         NewImagVal  =   absVal*sin(AngMatlab);
01347 
01348 
01349                         fkVecR[ kIy +kIx *End] =  NewRealVal ;
01350                         fkVecR[(End-1-kIy)+kCIx*End] =  NewRealVal ;
01351                         fkVecI[ kIy +kIx *End] =  NewImagVal ;
01352                         fkVecI[(End-1-kIy)+kCIx*End] = -NewImagVal ;
01353                         absD1fkVec[(End-1-kIy) + kIx  *End] = absVal;
01354                         absD1fkVec[(End-1-kIy) + kCIx *End] = absVal;
01355                         kVecX[kIy+kIx  *End] =  kx      ;
01356                         kVecX[kIy+kCIx *End] =  kCx    ;
01357                         kVecY[kIy+kIx  *End] =  ky     ;
01358                         kVecY[kIy+kCIx *End] =  ky     ;
01359 
01360  //                     cout << " kIxM= " << kIx+1 << " kIy=" << kIy+1 << " fkVecR[i] =" << NewRealVal << " fkVecI[i]="  << NewImagVal <<"  angle[i]= "  << AngMatlab << " Total Index" << kIy+kIx *End << endl;
01361 
01362 //                      printf("kx=%d,ky=%d,tempVal =%f+ i %4.2f \n",kx,ky,realVal,imagVal );
01363 //                      cout << "kx = " << kx << "; ky = "<< ky << "; val is" << realVal<<"+ i "<<imagVal<< endl;
01364 
01365 //                      cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; fkAng*9/ 2pi is " << fkAng*9/2/M_PI<<  endl;
01366 //                      cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; absval is " << absVal<<  "; realval is " << NewRealVal<< "; imagval is " << NewImagVal<< endl;
01367                 }
01368         }
01369 
01370 
01371         for (int TotalInd = 0 ;  TotalInd < CountxyMax ; TotalInd++){
01372                 int kx     = kVecX[TotalInd]; // This is the value of the index for a matlab image (-1)
01373                 int kIx    = kx+Mid-1; // This is the value of the index for a matlab image (-1)
01374                 int ky     = kVecY[TotalInd];
01375                 int kIy    = ky+Mid-1; // This is the value of the index for a matlab image (-1)
01376                 float fkR  = fkVecR[kIy+kIx *End]  ;
01377                 float fkI  = fkVecI[kIy+kIx *End]  ;
01378         }
01379 
01380         float frR= 3.0/4.0;
01381         frR= 1;
01382         int LradRange= (int) (1+floor(Mid/frR -.1)) ;
01383 
01384         float *radRange = new float[LradRange]; //= 0:.75:(Mid-1);
01385         for (int irad=0; irad < LradRange; irad++){
01386                         radRange[irad] =  frR*irad;
01387 //                      cout << " irad = " << irad << " radRange[irad]= " <<  radRange[irad] <<  " LradRange= " << LradRange << endl;
01388         }
01389 
01390         cout << "Starting the calculation of invariants" << endl;
01391 
01392 
01393         if (type==0) {
01394                 int LthetaRange  = 59;
01395                 float ftR        = (2.0f*M_PI/LthetaRange );
01396                 float *thetaRange = new float[LthetaRange]; //= 0:.75:(Mid-1);
01397 
01398                 for (int ith=0; ith < LthetaRange; ith++){
01399                                 thetaRange[ith] =  ftR*ith; }
01400 
01401 
01402 
01403 
01404 
01405                 int TotalVol = LradRange*LradRange*LthetaRange;
01406 
01407                 float *RotTransInv   = new  float[TotalVol];
01408                 float *WeightInv     = new  float[TotalVol];
01409 
01410                 for (int jW=0; jW<TotalVol; jW++) {
01411                         RotTransInv[jW] = 0;
01412                         WeightInv[jW]   = 0;
01413                 }
01414 
01415 
01416                 for (int jW=0; jW<TotalVol; jW++) {
01417                         RotTransInv[jW] = 0;
01418                         WeightInv[jW]   = 0;
01419                 }
01420         //      float  *RotTransInv       = new float[LradRange*LradRange ] ;
01421         //      float  *RotTransInvN      = new float[LradRange*LradRange*(NMax+1) ] ;
01422 
01423                 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){  // Main Section for type 0
01424                         int kx = kVecX[Countkxy] ;
01425                         int ky = kVecY[Countkxy] ;
01426                         float k2 = ::sqrt((float)(kx*kx+ky*ky));
01427                         float phiK =0;
01428                         if (k2>0)    phiK = jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1]; //  phiK=atan2(ky,kx);
01429                         float fkR     = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
01430                         float fkI     = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End]  ;
01431         //              printf("Countkxy=%d,\t kx=%d, ky=%d, fkR=%3.2f,fkI=%3.2f \n", Countkxy, kx, ky, fkR, fkI);
01432 
01433                         if ((k2==0)|| (k2>Mid) ) { continue;}
01434 
01435                         for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){   // This is the innermost loop
01436                                 int qx   = kVecX[Countqxy] ;
01437                                 int qy   = kVecY[Countqxy] ;
01438                                 float q2   = ::sqrt((float)(qx*qx+qy*qy));
01439                                 if ((q2==0)|| (q2>Mid) ) {continue;}
01440                                 float phiQ =0;
01441                                 if (q2>0)   phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1]; // phiQ=atan2(qy,qx);
01442                                 float fqR     = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
01443                                 float fqI     = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End]  ;
01444                                 int kCx  = (-kx-qx);
01445                                 int kCy  = (-ky-qy);
01446                                 int kCIx = ((kCx+Mid+2*End)%End);// labels of the image in C
01447                                 int kCIy = ((kCy+Mid+2*End)%End);
01448                                 kCx  = ((kCIx+End-1)%End)+1-Mid; // correct
01449                                 kCy  = ((kCIy+End-1)%End)+1-Mid ; // correct
01450 
01451                                 float C2   = ::sqrt((float)(kCx*kCx+ kCy*kCy));
01452                                 int CountCxy  = (kCx+Mid-1)*End+(kCy+Mid-1);
01453                                 float fCR     = fkVecR[CountCxy];
01454                                 float fCI     = fkVecI[CountCxy];
01455         /*                      if (Countkxy==1) {
01456                                         printf(" Countqxy=%d, absD1fkVec(Countqxy)=%f,qx=%d, qy=%d \n", Countqxy, absD1fkVec[Countqxy],qx, qy);
01457                                         printf(" CountCxy=%d, absD1fkVec[CountCxy]=%f,kCx=%d,kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
01458                                 }*/
01459                                 float   phiC = jxjyatan2[ (kCy+Mid-1)*End + kCx+Mid-1];
01460                                 float   phiQK = (4*M_PI+phiQ-phiK);
01461                                 while (phiQK> (2*M_PI)) phiQK -= (2*M_PI);
01462 
01463 
01464 
01465                                 float bispectemp  = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR  +fqR*fCI));
01466 
01467                                 if  ((q2<k2) )  continue;
01468 //                              if  ((q2<k2) || (C2<k2) || (C2<q2))  continue;
01469 
01470         //                              printf(" CountCxy=%d, absD1fkVec[CountCxy]=%f,kCx=%d,kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
01471 
01472         //                      up to here, matched perfectly with Matlab
01473 
01474                                 int k2IndLo  = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
01475                                 int k2IndHi = k2IndLo;
01476                                 float k2Lo= radRange[k2IndLo];
01477                                 if (k2IndLo+1< LradRange) {
01478                                         k2IndHi   = k2IndLo+1;
01479                                 }
01480                                 float k2Hi= radRange[k2IndHi];
01481 
01482                                 float kCof =k2-k2Lo;
01483 
01484                                 int q2IndLo  = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
01485                                 int q2IndHi=q2IndLo;
01486                                 float q2Lo= radRange[q2IndLo];
01487                                 if (q2IndLo+1 < LradRange)  {
01488                                         q2IndHi   = q2IndLo+1 ;
01489                                 }
01490                                 float qCof = q2-q2Lo;
01491 
01492                                 if ((qCof<0) || (qCof >1) ) {
01493                                         cout<< "Weird! qCof="<< qCof <<  " q2="<< q2 << " q2IndLo="<< q2IndLo << endl ;
01494                                         int x    ;
01495                                         cin >> x ;
01496                                 }
01497 
01498                                 int thetaIndLo = 0; while ((phiQK>=thetaRange[thetaIndLo+1])&& (thetaIndLo+1<LthetaRange)) thetaIndLo +=1;
01499                                 int thetaIndHi = thetaIndLo;
01500 
01501                                 float thetaLo  = thetaRange[thetaIndLo];
01502                                 float thetaHi = thetaLo;
01503                                 float thetaCof = 0;
01504 
01505                                 if (thetaIndLo+1< LthetaRange) {
01506                                         thetaIndHi = thetaIndLo +1;
01507                                 }else{
01508                                         thetaIndHi=0;
01509                                 }
01510 
01511                                 thetaHi    = thetaRange[thetaIndHi];
01512 
01513                                 if (thetaHi==thetaLo) {
01514                                         thetaCof =0 ;
01515                                 } else {
01516                                         thetaCof   = (phiQK-thetaLo)/(thetaHi-thetaLo);
01517                                 }
01518 
01519                                 if ((thetaCof>2*M_PI)  ) {
01520                                         cout<< "Weird! thetaCof="<< thetaCof <<endl ;
01521                                         thetaCof=0;
01522                                 }
01523 
01524 
01525         //                      if ((thetaIndLo>=58) || (k2IndLo >= LradRange-1) || (q2IndLo >= LradRange-1) ) {
01526 
01527 
01528                                 for (int jk =1; jk<=2; jk++){
01529                                 for (int jq =1; jq<=2; jq++){
01530                                 for (int jtheta =1; jtheta<=2; jtheta++){
01531 
01532                                         float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1))
01533                                                         * (thetaCof+(1-2*thetaCof)*(jtheta==1));
01534 
01535 
01536                                         int k2Ind      =  k2IndLo*(jk==1)      +   k2IndHi*(jk==2);
01537                                         int q2Ind      =  q2IndLo*(jq==1)      +   q2IndHi*(jq==2);
01538                                         int thetaInd   =  thetaIndLo*(jtheta==1)  + thetaIndHi*(jtheta ==2);
01539                                         int TotalInd   = thetaInd*LradRange*LradRange+q2Ind*LradRange+k2Ind;
01540         /*                              if (TotalInd+1 >=  LthetaRange*LradRange*LradRange) {
01541                                                 cout << "Weird!!! TotalInd="<< TotalInd << " IndMax" << LthetaRange*LradRange*LradRange << " LradRange=" << LradRange << endl;
01542                                                 cout << "k2Ind= "<< k2Ind  << " q2Ind="<< q2Ind  << " thetaInd="<< thetaInd  << " q2IndLo="<< q2IndLo  << " q2IndHi="<< q2IndHi  <<  endl;
01543                                                 cout << "k2=" << k2 << "q2=" << q2 << " phiQK=" << phiQK*180.0/M_PI<< endl;
01544                                         }*/
01545 
01546                                         RotTransInv[TotalInd] += Weight*bispectemp;
01547                                         WeightInv[TotalInd]   +=  Weight;
01548         //                              cout << "k2Ind= "<< k2Ind  << " q2Ind="<< q2Ind  << "Weight=" << Weight << endl;
01549                                 }}}
01550                         } // Countqxy
01551                 } // Countkxy
01552 
01553                 cout << "Finished Main Section " << endl;
01554 
01555         /*              RotTransInvN[jr1 + LradRange*jr2+LradRange*LradRange*N] = RotTransInvTemp  ;*/
01556 
01557                 cout << " LradRange " <<LradRange <<" LthetaRange " << LthetaRange << endl;
01558                 EMData *RotTransInvF  = new  EMData(LradRange,LradRange,LthetaRange);
01559                 EMData *WeightImage   = new  EMData(LradRange,LradRange,LthetaRange);
01560 
01561         //      cout << "FFFFFFF" << endl;
01562         //
01563         //      RotTransInvF -> set_size(LradRange,LradRange,LthetaRange);
01564         //
01565         //      cout << "GGGG" << endl;
01566 
01567                 for (int jtheta =0; jtheta < LthetaRange; jtheta++){    // write out to RotTransInvF
01568                 for (int jq =0; jq<LradRange; jq++){ // LradRange
01569                 for (int jk =0; jk<LradRange ; jk++){// LradRange
01570         //              cout << "Hi There" << endl;
01571                         int TotalInd   = jtheta*LradRange*LradRange+jq*LradRange+jk;
01572                         float Weight = WeightInv[TotalInd];
01573                         WeightImage    -> set_value_at(jk,jq,jtheta,Weight);
01574                         RotTransInvF   -> set_value_at(jk,jq,jtheta,0);
01575                         if (Weight <= 0) continue;
01576                         RotTransInvF -> set_value_at(jk,jq,jtheta,RotTransInv[TotalInd] / Weight);//  include /Weight
01577                         int newjtheta = (LthetaRange- jtheta)%LthetaRange;
01578                         RotTransInvF -> set_value_at(jq,jk,newjtheta,RotTransInv[TotalInd]/Weight );//  include /Weight
01579                                 }
01580                         }
01581                 }
01582 
01583                 cout << " Almost Done " << endl;
01584                 system("rm -f WeightImage.???");
01585                 WeightImage  -> write_image("WeightImage.img");
01586 
01587                 return  RotTransInvF ;
01588         } // Finish type 0
01589 
01590         if (type==1) {
01591                 int TotalVol = LradRange*LradRange;
01592 
01593                 float *RotTransInv   = new  float[TotalVol];
01594                 float *WeightInv     = new  float[TotalVol];
01595 
01596                 for (int jW=0; jW<TotalVol; jW++) {
01597                         RotTransInv[jW] = 0;
01598                         WeightInv[jW]   = 0;
01599                 }
01600 
01601 
01602                 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){
01603                         int kx = kVecX[Countkxy] ;
01604                         int ky = kVecY[Countkxy] ;
01605                         float k2 = ::sqrt((float)(kx*kx+ky*ky));
01606                         float fkR     = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
01607                         float fkI     = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End]  ;
01608         //              printf("Countkxy=%d,\t kx=%d, ky=%d, fkR=%3.2f,fkI=%3.2f \n", Countkxy, kx, ky, fkR, fkI);
01609 
01610                         if ((k2==0)|| (k2>Mid) ) { continue;}
01611 
01612                         for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){   // This is the innermost loop
01613 
01614 //                      up to here, matched perfectly with Matlab
01615                                 int qx   = kVecX[Countqxy] ;
01616                                 int qy   = kVecY[Countqxy] ;
01617                                 float q2   = ::sqrt((float)(qx*qx+qy*qy));
01618                                 if ((q2==0)|| (q2>Mid) ) {continue;}
01619                                 if  ((q2<k2) )   continue;
01620 
01621                                 float fqR     = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
01622                                 float fqI     = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End]  ;
01623 
01624                                 int kCx  = (-kx-qx);
01625                                 int kCy  = (-ky-qy);
01626                                 int kCIx = ((kCx+Mid+2*End)%End);// labels of the image in C
01627                                 int kCIy = ((kCy+Mid+2*End)%End);
01628                                 kCx  = ((kCIx+End-1)%End)+1-Mid; // correct
01629                                 kCy  = ((kCIy+End-1)%End)+1-Mid ; // correct
01630 
01631                                 float C2   = ::sqrt((float)(kCx*kCx+ kCy*kCy));
01632                                 int CountCxy  = (kCx+Mid-1)*End+(kCy+Mid-1);
01633                                 float fCR     = fkVecR[CountCxy];
01634                                 float fCI     = fkVecI[CountCxy];
01635 
01636 
01637                                 float bispectemp  = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR  +fqR*fCI));
01638 
01639 
01640                                 int k2IndLo  = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
01641                                 int k2IndHi = k2IndLo;
01642                                 float k2Lo= radRange[k2IndLo];
01643                                 if (k2IndLo+1< LradRange) {
01644                                         k2IndHi   = k2IndLo+1;
01645                                 }
01646                                 float k2Hi= radRange[k2IndHi];
01647 
01648                                 float kCof =k2-k2Lo;
01649 
01650                                 int q2IndLo  = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
01651                                 int q2IndHi=q2IndLo;
01652                                 float q2Lo= radRange[q2IndLo];
01653                                 if (q2IndLo+1 < LradRange)  {
01654                                         q2IndHi   = q2IndLo+1 ;
01655                                 }
01656                                 float qCof = q2-q2Lo;
01657 
01658 
01659                                 for (int jk =1; jk<=2; jk++){
01660                                 for (int jq =1; jq<=2; jq++){
01661 
01662                                         float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1));
01663 
01664                                         int k2Ind      =  k2IndLo*(jk==1)      +   k2IndHi*(jk==2);
01665                                         int q2Ind      =  q2IndLo*(jq==1)      +   q2IndHi*(jq==2);
01666                                         int TotalInd   = q2Ind*LradRange+k2Ind;
01667                                         RotTransInv[TotalInd] += Weight*bispectemp;
01668                                         WeightInv[TotalInd]   +=  Weight;
01669         //                              cout << "k2Ind= "<< k2Ind  << " q2Ind="<< q2Ind  << "Weight=" << Weight << endl;
01670                                 }}
01671                         } // Countqxy
01672                 } // Countkxy
01673 
01674 //              cout << "Finished Main Section " << endl;
01675 //              cout << " LradRange " <<LradRange <<  endl;
01676 
01677 
01678                 EMData *RotTransInvF  = new  EMData(LradRange,LradRange);
01679                 EMData *WeightImage   = new  EMData(LradRange,LradRange);
01680 
01681                 for (int jk =0; jk<LradRange ; jk++){// LradRange
01682                 for (int jq =jk; jq<LradRange; jq++){ // LradRange
01683                         int TotalInd      = jq*LradRange+jk;
01684                         int TotalIndBar   = jq*LradRange+jk;
01685                         float Weight = WeightInv[TotalInd] + WeightInv[TotalIndBar];
01686                         if (Weight <=0) continue;
01687                         WeightImage    -> set_value_at(jk,jq,Weight);
01688                         WeightImage    -> set_value_at(jq,jk,Weight);
01689 #ifdef _WIN32
01690                         float ValNow  = pow( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight, 1.0f/3.0f )  ;
01691 #else
01692                         float ValNow  = cbrt( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight )  ;
01693 #endif  //_WIN32
01694                         RotTransInvF -> set_value_at(jk,jq,ValNow);//  include /Weight
01695                         RotTransInvF -> set_value_at(jq,jk,ValNow );//  include /Weight
01696                 }}
01697 
01698                 system("rm -f WeightImage.???");
01699                 WeightImage  -> write_image("WeightImage.img");
01700 
01701                 return  RotTransInvF ;
01702         }
01703 }
01704 
01705 
01706 void EMData::insert_clip(const EMData * const block, const IntPoint &origin) {
01707         int nx1 = block->get_xsize();
01708         int ny1 = block->get_ysize();
01709         int nz1 = block->get_zsize();
01710 
01711         Region area(origin[0], origin[1], origin[2],nx1, ny1, nz1);
01712 
01713         int x0 = (int) area.origin[0];
01714         x0 = x0 < 0 ? 0 : x0;
01715 
01716         int y0 = (int) area.origin[1];
01717         y0 = y0 < 0 ? 0 : y0;
01718 
01719         int z0 = (int) area.origin[2];
01720         z0 = z0 < 0 ? 0 : z0;
01721 
01722         int x1 = (int) (area.origin[0] + area.size[0]);
01723         x1 = x1 > nx ? nx : x1;
01724 
01725         int y1 = (int) (area.origin[1] + area.size[1]);
01726         y1 = y1 > ny ? ny : y1;
01727 
01728         int z1 = (int) (area.origin[2] + area.size[2]);
01729         z1 = z1 > nz ? nz : z1;
01730         if (z1 <= 0) {
01731                 z1 = 1;
01732         }
01733 
01734         int xd0 = (int) (area.origin[0] < 0 ? -area.origin[0] : 0);
01735         int yd0 = (int) (area.origin[1] < 0 ? -area.origin[1] : 0);
01736         int zd0 = (int) (area.origin[2] < 0 ? -area.origin[2] : 0);
01737 
01738         if (x1 < x0 || y1 < y0 || z1 < z0) return; // out of bounds, this is fine, nothing happens
01739 
01740         size_t clipped_row_size = (x1-x0) * sizeof(float);
01741         int src_secsize =  nx1 * ny1;
01742         int dst_secsize = nx * ny;
01743 
01744 #ifdef EMAN2_USING_CUDA
01745         if (gpu_operation_preferred()) {
01746 //              cudaMemcpy3DParms copyParams = {0};
01747 //              copyParams.srcPtr   = make_cudaPitchedPtr((void*)get_cuda_data(), nx*sizeof(float), nx, ny);
01748 //              cout << "Src pos is " << x0 << " " << y0 << " " << z0 << endl;
01749 //              copyParams.srcPos  = make_cudaPos(x0,y0,z0);
01750 //              cout << "Extent is " << x1-x0 << " " << y1-y0 << " " << z1-z0 << endl;
01751 //              cudaExtent extent = make_cudaExtent(x1-x0,y1-y0,z1-z0);
01752 //              int rnx = result->get_xsize();
01753 //              int rny = result->get_ysize();
01754 //              copyParams.dstPtr = make_cudaPitchedPtr((void*)result->get_cuda_data(),rnx*sizeof(float),rnx,rny );
01755 //              cout << "Dest pos is " << xd0 << " " << yd0 << " " << zd0 << endl;
01756 //              copyParams.dstPos  = make_cudaPos(xd0,yd0,zd0);
01757 //              copyParams.extent   = extent;
01758 //              copyParams.kind     = cudaMemcpyDeviceToDevice;
01759 //              cudaError_t error =  cudaMemcpy3D(&copyParams);
01760 //              if ( error != cudaSuccess) {
01761 //                      string e = cudaGetErrorString(error);
01762 //                      throw UnexpectedBehaviorException("CudaMemcpy3D failed with error: " + e);
01763 //              }
01764                 get_cuda_data();
01765                 CudaDataLock lock(this);
01766                 block->get_cuda_data();
01767                 CudaDataLock lock2(block);
01768                 for (int i = 0; i < (z1-z0); i++) {
01769                         float* ldst = get_cuda_data() + (z0+i) * dst_secsize + y0 * nx + x0;
01770                         float* lsrc = block->get_cuda_data() + (zd0+i) * src_secsize + yd0 * nx1 + xd0;
01771 
01772 
01773 //                      float* ldst = result->get_cuda_data() + (zd0+i) * dst_secsize + yd0 * (int)area.size[0] + xd0;
01774 //                      float* lsrc = get_cuda_data() + (z0+i) * src_secsize + y0 * nx + x0;
01775                         cudaError_t error =  cudaMemcpy2D( ldst, (nx)*sizeof(float), lsrc,nx1*sizeof(float), clipped_row_size, (y1-y0), cudaMemcpyDeviceToDevice );
01776                         if ( error != cudaSuccess) {
01777                                 string e = cudaGetErrorString(error);
01778                                 throw UnexpectedBehaviorException("cudaMemcpy2D failed in get_clip with error: " + e);
01779                         }
01780                 }
01781                 gpu_update();
01782                 EXITFUNC;
01783                 return;
01784         }
01785 #endif
01786         float *src = block->get_data() + zd0 * src_secsize + yd0 * nx1 + xd0;
01787         float *dst = get_data() + z0 * dst_secsize + y0 * nx + x0;
01788 
01789 //              float *src = get_data() + z0 * src_secsize + y0 * nx + x0;
01790 //              float *dst = result->get_data();
01791 //              dst += zd0 * dst_secsize + yd0 * (int)area.size[0] + xd0;
01792 
01793         int src_gap = src_secsize - (y1-y0) * nx1;
01794         int dst_gap = dst_secsize - (y1-y0) * nx;
01795         for (int i = z0; i < z1; i++) {
01796                 for (int j = y0; j < y1; j++) {
01797                         EMUtil::em_memcpy(dst, src, clipped_row_size);
01798                         src += nx1;
01799                         dst += nx;
01800                 }
01801                 src += src_gap;
01802                 dst += dst_gap;
01803         }
01804 
01805         update();
01806         EXITFUNC;
01807 }
01808 
01809 
01810 void EMData::insert_scaled_sum(EMData *block, const FloatPoint &center,
01811                                                    float scale, float)
01812 {
01813         ENTERFUNC;
01814         float * data = get_data();
01815         if (get_ndim()==3) {
01816                 // Start by determining the region to operate on
01817                 int xs=(int)floor(block->get_xsize()*scale/2.0);
01818                 int ys=(int)floor(block->get_ysize()*scale/2.0);
01819                 int zs=(int)floor(block->get_zsize()*scale/2.0);
01820                 int x0=(int)center[0]-xs;
01821                 int x1=(int)center[0]+xs;
01822                 int y0=(int)center[1]-ys;
01823                 int y1=(int)center[1]+ys;
01824                 int z0=(int)center[2]-zs;
01825                 int z1=(int)center[2]+zs;
01826 
01827                 if (x1<0||y1<0||z1<0||x0>get_xsize()||y0>get_ysize()||z0>get_zsize()) return;   // object is completely outside the target volume
01828 
01829                 // make sure we stay inside the volume
01830                 if (x0<0) x0=0;
01831                 if (y0<0) y0=0;
01832                 if (z0<0) z0=0;
01833                 if (x1>=get_xsize()) x1=get_xsize()-1;
01834                 if (y1>=get_ysize()) y1=get_ysize()-1;
01835                 if (z1>=get_zsize()) z1=get_zsize()-1;
01836 
01837                 float bx=block->get_xsize()/2.0f;
01838                 float by=block->get_ysize()/2.0f;
01839                 float bz=block->get_zsize()/2.0f;
01840 
01841                 size_t idx;
01842                 for (int x=x0; x<=x1; x++) {
01843                         for (int y=y0; y<=y1; y++) {
01844                                 for (int z=z0; z<=z1; z++) {
01845                                         idx = x + y * nx + z * nx * ny;
01846                                         data[idx] +=
01847                                                 block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by,(z-center[2])/scale+bz);
01848                                 }
01849                         }
01850                 }
01851                 update();
01852         }
01853         else if (get_ndim()==2) {
01854                 // Start by determining the region to operate on
01855                 int xs=(int)floor(block->get_xsize()*scale/2.0);
01856                 int ys=(int)floor(block->get_ysize()*scale/2.0);
01857                 int x0=(int)center[0]-xs;
01858                 int x1=(int)center[0]+xs;
01859                 int y0=(int)center[1]-ys;
01860                 int y1=(int)center[1]+ys;
01861 
01862                 if (x1<0||y1<0||x0>get_xsize()||y0>get_ysize()) return; // object is completely outside the target volume
01863 
01864                 // make sure we stay inside the volume
01865                 if (x0<0) x0=0;
01866                 if (y0<0) y0=0;
01867                 if (x1>=get_xsize()) x1=get_xsize()-1;
01868                 if (y1>=get_ysize()) y1=get_ysize()-1;
01869 
01870                 float bx=block->get_xsize()/2.0f;
01871                 float by=block->get_ysize()/2.0f;
01872 
01873                 for (int x=x0; x<=x1; x++) {
01874                         for (int y=y0; y<=y1; y++) {
01875                                 data[x + y * nx] +=
01876                                         block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by);
01877                         }
01878                 }
01879                 update();
01880         }
01881         else {
01882                 LOGERR("insert_scaled_sum supports only 2D and 3D data");
01883                 throw ImageDimensionException("2D/3D only");
01884         }
01885 
01886         EXITFUNC;
01887 }
01888 //                      else if ( m == 0 )
01889 //                      {
01890 //                              if ( n_f == -ny/2 )
01891 //                              {
01892 //                                      t2++;
01893 // //                                   continue;
01894 //                                      for (int y = 0; y < return_slice->get_ysize(); ++y) {
01895 //                                              for (int x = 0; x < return_slice->get_xsize(); ++x) {
01896 //                                                      double cur_val = return_slice->get_value_at(x,y);
01897 //                                                      return_slice->set_value_at(x,y,cur_val+dat[idx]*std::pow(-1.0f,y));
01898 //                                              }
01899 //                                      }
01900 //                                      if (phase > 0.01 ) cout << "foo 2 " << phase << " " << amp << " " << dat[idx] << endl;
01901 //                              }
01902 //                              else
01903 //                              {
01904 //                                      if ( n_f < 1 ) continue;
01905 //                                      t3++;
01906 //                                      for (int y = 0; y < return_slice->get_ysize(); ++y) {
01907 //                                              for (int x = 0; x < return_slice->get_xsize(); ++x) {
01908 //                                                      double cur_val = return_slice->get_value_at(x,y);
01909 //                                                      return_slice->set_value_at(x,y,cur_val+2*amp*cos(ndash*y+phase));
01910 //                                              }
01911 //                                      }
01912 //                              }
01913 //                      }
01914 //                      else if ( n_f == -ny/2 )
01915 //                      {
01916 //                              if ( m == ((nx-2)/2) )
01917 //                              {
01918 //                                      t4++;
01919 //                                      for (int y = 0; y < return_slice->get_ysize(); ++y) {
01920 //                                              for (int x = 0; x < return_slice->get_xsize(); ++x) {
01921 //                                                      double cur_val = return_slice->get_value_at(x,y);
01922 //                                                      return_slice->set_value_at(x,y,cur_val+dat[idx]*std::pow(-1.0f,x+y));
01923 //                                              }
01924 //                                      }
01925 //                                      if (phase > 0.01 ) cout << "foo 4 " << phase << " " << amp << " " << dat[idx] << endl;
01926 //                              }
01927 //                              else
01928 //                              {
01929 //                                      t5++;
01930 //                                      for (int y = 0; y < return_slice->get_ysize(); ++y) {
01931 //                                              for (int x = 0; x < return_slice->get_xsize(); ++x) {
01932 //                                                      double cur_val = return_slice->get_value_at(x,y);
01933 //                                                      return_slice->set_value_at(x,y,cur_val+2*amp*cos(mdash*x+phase));
01934 //                                              }
01935 //                                      }
01936 //                              }
01937 //                      }
01938 //                      else if ( n_f == 0 )
01939 //                      {
01940 //                              if ( m == ((nx-2)/2) )
01941 //                              {
01942 //                                      t6++;
01943 //                                      for (int y = 0; y < return_slice->get_ysize(); ++y) {
01944 //                                              for (int x = 0; x < return_slice->get_xsize(); ++x) {
01945 //                                                      double cur_val = return_slice->get_value_at(x,y);
01946 //                                                      return_slice->set_value_at(x,y,cur_val+dat[idx]*std::pow(-1.0f,x));
01947 //                                              }
01948 //                                      }
01949 //                                      if (phase > 0.01 ) cout << "foo 3 " << phase << " " << amp << " " << dat[idx] << endl;
01950 //                              }
01951 //                              else
01952 //                              {
01953 //                                      t7++;
01954 //                                      for (int y = 0; y < return_slice->get_ysize(); ++y) {
01955 //                                              for (int x = 0; x < return_slice->get_xsize(); ++x) {
01956 //                                                      double cur_val = return_slice->get_value_at(x,y);
01957 //                                                      return_slice->set_value_at(x,y,cur_val+2*amp*cos(mdash*x+M_PI*y + phase));
01958 //                                              }
01959 //                                      }
01960 //                              }
01961 //                      }
01962 //                      else if ( m == ((nx-2)/2) )
01963 //                      {
01964 //                              if ( n_f < 1 ) continue;
01965 //                              t8++;
01966 //                              for (int y = 0; y < return_slice->get_ysize(); ++y) {
01967 //                                      for (int x = 0; x < return_slice->get_xsize(); ++x) {
01968 //                                              double cur_val = return_slice->get_value_at(x,y);
01969 //                                              return_slice->set_value_at(x,y,cur_val+2*amp*cos(ndash*y+M_PI*x+phase));
01970 //                                      }
01971 //                              }
01972 //                      }
01973 // }

Generated on Fri Apr 30 15:38:51 2010 for EMAN2 by  doxygen 1.3.9.1