EMAN2
emdata_transform.cpp
Go to the documentation of this file.
1/*
2 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
3 * Copyright (c) 2000-2006 Baylor College of Medicine
4 *
5 * This software is issued under a joint BSD/GNU license. You may use the
6 * source code in this file under either license. However, note that the
7 * complete EMAN2 and SPARX software packages have some GPL dependencies,
8 * so you are responsible for compliance with the licenses of these packages
9 * if you opt to use BSD licensing. The warranty disclaimer below holds
10 * in either instance.
11 *
12 * This complete copyright notice must be included in any revised version of the
13 * source code. Additional authorship citations may be added, but existing
14 * author citations must be preserved.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 *
30 * */
31
32#include "emdata.h"
33#include "emfft.h"
34
35#include <cstring>
36#include <cstdio>
37
38#include "gsl/gsl_sf_result.h"
39#include "gsl/gsl_sf_bessel.h"
40#include <iostream>
41#include <algorithm>
42#include <vector>
43#include <utility>
44#include <cmath>
45#include "util.h"
46
47//#ifdef EMAN2_USING_CUDA
48//#include "cuda/cuda_processor.h"
49//#endif
50
51using namespace EMAN;
52using namespace std;
53typedef vector< pair<float,int> > vp;
54
55EMData *EMData::do_fft() const
56{
57 ENTERFUNC;
58#ifdef FFT_CACHING
59 if (fftcache!=0) {
60 return fftcache->copy();
61 }
62#endif //FFT_CACHING
63
64 if (is_complex() ) { // ming add 08/17/2010
65#ifdef NATIVE_FFT
66 LOGERR(" NATIVE_FFT does not support complex to complex."); // PAP
67 throw ImageFormatException("real image expected. Input image is complex image.");
68#else
69 EMData *temp_in=copy();
70 EMData *dat= copy_head();
71 int offset;
72 if(is_fftpadded()) {
73 offset = is_fftodd() ? 1 : 2;
74 }
75 else offset=0;
76 //printf("offset=%d\n",offset);
77 EMfft::complex_to_complex_nd(temp_in->get_data(),dat->get_data(),nx-offset,ny,nz);
78
79 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(true);
80
81 dat->update();
82 delete temp_in;
83 EXITFUNC;
84 return dat;
85#endif // NATIVE_FFT
86 } else {
87 int nxreal = nx;
88 int offset = 2 - nx%2;
89 int nx2 = nx + offset;
90 EMData* dat = copy_head();
91 dat->set_size(nx2, ny, nz);
92 //dat->to_zero(); // do not need it, real_to_complex will do it right anyway
93 if (offset == 1) dat->set_fftodd(true);
94 else dat->set_fftodd(false);
95
96 float *d = dat->get_data();
97 //std::cout<<" do_fft "<<rdata[5]<<" "<<d[5]<<std::endl;
98 EMfft::real_to_complex_nd(get_data(), d, nxreal, ny, nz);
99
100 dat->update();
101 dat->set_fftpad(true);
102 dat->set_complex(true);
103 dat->set_attr("is_intensity",false);
104 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(true);
105 dat->set_ri(true);
106
107 EXITFUNC;
108#ifdef FFT_CACHING
109// printf("%p %d\n",this,nxyz);
110 if (nxyz<80000000) {
111 fftcache=dat->copy();
112 }
113#endif //FFT_CACHING
114 return dat;
115 }
116}
117
118void EMData::do_fft_inplace()
119{
120 ENTERFUNC;
121
122 if ( is_complex() ) {
123 LOGERR("real image expected. Input image is complex image.");
124 throw ImageFormatException("real image expected. Input image is complex image.");
125 }
126
127 size_t offset;
128 int nxreal;
129 get_data(); // Required call if GPU caching is being used. Otherwise harmless
130 if (!is_fftpadded()) {
131 // need to extend the matrix along x
132 // meaning nx is the un-fftpadded size
133 nxreal = nx;
134 offset = 2 - nx%2;
135 if (1 == offset) set_fftodd(true);
136 else set_fftodd(false);
137 int nxnew = nx + offset;
138 set_size(nxnew, ny, nz);
139 for (int iz = nz-1; iz >= 0; iz--) {
140 for (int iy = ny-1; iy >= 0; iy--) {
141 for (int ix = nxreal-1; ix >= 0; ix--) {
142 size_t oldxpos = ix + (iy + iz*ny)*(size_t)nxreal;
143 size_t newxpos = ix + (iy + iz*ny)*(size_t)nxnew;
144 (*this)(newxpos) = (*this)(oldxpos);
145 }
146 }
147 }
148 set_fftpad(true);
149 } else {
150 offset = is_fftodd() ? 1 : 2;
151 nxreal = nx - offset;
152 }
153 EMfft::real_to_complex_nd(rdata, rdata, nxreal, ny, nz);
154
155 set_complex(true);
156 if(ny==1 && nz==1) set_complex_x(true);
157 set_ri(true);
158
159 update();
160
161 EXITFUNC;
162 //return this;
163}
164
165#ifdef EMAN2_USING_CUDA
166
167#include "cuda/cuda_emfft.h"
168
169EMData *EMData::do_fft_cuda()
170{
171 ENTERFUNC;
172
173 if ( is_complex() ) {
174 LOGERR("real image expected. Input image is complex image.");
175 throw ImageFormatException("real image expected. Input image is complex image.");
176 }
177
178 int offset = 2 - nx%2;
179 EMData* dat = new EMData(0,0,nx+offset,ny,nz,attr_dict);
180 if(!dat->rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");
181 //cout << "Doing CUDA FFT " << cudarwdata << endl;
182 if(cudarwdata == 0){copy_to_cuda();}
183 cuda_dd_fft_real_to_complex_nd(cudarwdata, dat->cudarwdata, nx, ny,nz, 1);
184
185 if (offset == 1) dat->set_fftodd(true);
186 else dat->set_fftodd(false);
187
188 dat->set_fftpad(true);
189 dat->set_complex(true);
190 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(true);
191 dat->set_ri(true);
192 dat->update();
193
194 EXITFUNC;
195 return dat;
196}
197
198void EMData::do_fft_inplace_cuda()
199{
200 ENTERFUNC;
201
202 if ( is_complex() ) {
203 LOGERR("real image expected. Input image is complex image.");
204 throw ImageFormatException("real image expected. Input image is complex image.");
205 }
206
207 int offset = 2 - nx%2;
208 float* tempcudadata = 0;
209 cudaError_t error = cudaMalloc((void**)&tempcudadata,(nx + offset)*ny*nz*sizeof(float));
210 if( error != cudaSuccess) throw ImageFormatException("Couldn't allocate memory.");
211
212 //cout << "Doing CUDA FFT inplace" << cudarwdata << endl;
213 if(cudarwdata == 0){copy_to_cuda();}
214 cuda_dd_fft_real_to_complex_nd(cudarwdata, tempcudadata, nx, ny,nz, 1);
215 // this section is a bit slight of hand it actually does the FFT out of place but this avoids and EMData object creation and detruction...
216 cudaError_t ferror = cudaFree(cudarwdata);
217 if ( ferror != cudaSuccess) throw UnexpectedBehaviorException( "CudaFree failed:" + string(cudaGetErrorString(error)));
218 cudarwdata = tempcudadata;
219 num_bytes = (nx + offset)*ny*nz*sizeof(float);
220
221 if (offset == 1) set_fftodd(true);
222 else set_fftodd(false);
223
224 nx = nx + offset; // don't want to call set_size b/c that will delete my cudadata, remember what I am doing is a bit slignt of hand....
225 set_fftpad(true);
226 set_complex(true);
227 if(get_ysize()==1 && get_zsize()==1) set_complex_x(true);
228 set_ri(true);
229 update();
230
231 EXITFUNC;
232// return this;
233}
234
235EMData *EMData::do_ift_cuda()
236{
237 ENTERFUNC;
238
239 if (!is_complex()) {
240 LOGERR("complex image expected. Input image is real image.");
241 throw ImageFormatException("complex image expected. Input image is real image.");
242 }
243
244 if (!is_ri()) {
245 throw ImageFormatException("complex ri expected. Got amplitude/phase.");
246 }
247
248 int offset = is_fftodd() ? 1 : 2;
249 EMData* dat = new EMData(0,0,nx-offset,ny,nz,attr_dict);
250 if(!dat->rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");
251
252 if(cudarwdata == 0){copy_to_cuda();}
253
254
255 int ndim = get_ndim();
256 if ( ndim == 1 ) {
257 cuda_dd_fft_complex_to_real_nd(cudarwdata,dat->cudarwdata, nx-offset,1,1,1);
258 } else if (ndim == 2) {
259 cuda_dd_fft_complex_to_real_nd(cudarwdata,dat->cudarwdata, ny,nx-offset,1,1);
260 } else if (ndim == 3) {
261 cuda_dd_fft_complex_to_real_nd(cudarwdata,dat->cudarwdata, nz,ny,nx-offset,1);
262 } else throw ImageDimensionException("No cuda FFT support of images with dimensions exceeding 3");
263
264 // SCALE the inverse FFT
265 float scale = 1.0f/static_cast<float>((dat->get_size()));
266 dat->mult(scale);
267
268 dat->set_fftpad(false);
269 dat->set_fftodd(false);
270 dat->set_complex(false);
271 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(false);
272 dat->set_ri(false);
273// dat->gpu_update();
274 dat->update();
275
276 EXITFUNC;
277 return dat;
278}
279
280/*
281 FFT in place does not depad, hence this routine is of limited use b/c mem operations on the device are quite SLOW, JFF
282 use
283*/
284
285void EMData::do_ift_inplace_cuda()
286{
287 ENTERFUNC;
288
289 if (!is_complex()) {
290 LOGERR("complex image expected. Input image is real image.");
291 throw ImageFormatException("complex image expected. Input image is real image.");
292 }
293
294 if (!is_ri()) {
295 LOGWARN("run IFT on AP data, only RI should be used. ");
296 }
297
298 int offset = is_fftodd() ? 1 : 2;
299
300 if(cudarwdata == 0){copy_to_cuda();}
301
302 int ndim = get_ndim();
303 if ( ndim == 1 ) {
304 cuda_dd_fft_complex_to_real_nd(cudarwdata,cudarwdata, nx-offset,1,1,1);
305 } else if (ndim == 2) {
306 cuda_dd_fft_complex_to_real_nd(cudarwdata,cudarwdata, ny,nx-offset,1,1);
307 } else if (ndim == 3) {
308 cuda_dd_fft_complex_to_real_nd(cudarwdata,cudarwdata, nz,ny,nx-offset,1);
309 } else throw ImageDimensionException("No cuda FFT support of images with dimensions exceeding 3");
310#if defined USE_FFTW3 //native fft and ACML already done normalization
311 // SCALE the inverse FFT
312 int nxo = nx - offset;
313 float scale = 1.0f / (nxo * ny * nz);
314 mult(scale); //if we are just doing a CCF, this is a waste!
315#endif //USE_FFTW3
316
317 set_fftpad(true);
318 set_complex(false);
319
320 if(ny==1 && nz==1) set_complex_x(false);
321 set_ri(false);
322 update();
323
324 EXITFUNC;
325// return this;
326}
327
328#endif //EMAN2_USING_CUDA
329
330EMData *EMData::do_ift()
331{
332 ENTERFUNC;
333
334 if (!is_complex()) {
335 LOGERR("complex image expected. Input image is real image.");
336 throw ImageFormatException("complex image expected. Input image is real image.");
337 }
338
339 if (!is_ri()) {
340 LOGWARN("run IFT on AP data, only RI should be used. Converting.");
341 }
342
343 get_data(); // Required call if GPU caching is being used. Otherwise harmless
344 EMData* dat = copy_head();
345 dat->set_size(nx, ny, nz);
346 ap2ri();
347
348 float *d = dat->get_data();
349 int ndim = get_ndim();
350
351 /* Do inplace IFT on a image copy, because the complex to real transform of
352 * nd will destroy its input array even for out-of-place transforms.
353 */
354 memcpy((char *) d, (char *) rdata, (size_t)nx * ny * nz * sizeof(float));
355
356 int offset = is_fftodd() ? 1 : 2;
357 //cout << "Sending offset " << offset << " " << nx-offset << endl;
358 if (ndim == 1) {
359 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
360 } else {
361 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
362
363 size_t row_size = (nx - offset) * sizeof(float);
364 for (size_t i = 1; i < (size_t)ny * nz; i++) {
365 memmove((char *) &d[i * (nx - offset)], (char *) &d[i * nx], row_size);
366 }
367 }
368
369 dat->set_size(nx - offset, ny, nz); //remove the padding
370#if defined USE_FFTW3 //native fft and ACML already done normalization
371 // SCALE the inverse FFT
372 float scale = 1.0f / ((nx - offset) * ny * nz);
373 dat->mult(scale);
374#endif //USE_FFTW3
375 dat->set_fftodd(false);
376 dat->set_fftpad(false);
377 dat->set_complex(false);
378 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(false);
379 dat->set_ri(false);
380 dat->set_attr("is_intensity",false);
381 dat->update();
382
383
384 EXITFUNC;
385 return dat;
386}
387
388/*
389 FFT in place does not depad, return real x-extended image (needs to be depadded before use as PAP does in CCF routines)
390 use
391*/
392void EMData::do_ift_inplace()
393{
394 ENTERFUNC;
395
396 if (!is_complex()) {
397 LOGERR("complex image expected. Input image is real image.");
398 throw ImageFormatException("complex image expected. Input image is real image.");
399 }
400
401 if (!is_ri()) {
402 LOGWARN("run IFT on AP data, only RI should be used. ");
403 }
404 ap2ri();
405
406 int offset = is_fftodd() ? 1 : 2;
407 float* data = get_data();
408 EMfft::complex_to_real_nd(data, data, nx - offset, ny, nz);
409
410#if defined USE_FFTW3 //native fft and ACML already done normalization
411 // SCALE the inverse FFT
412 int nxo = nx - offset;
413 float scale = 1.0f / ((size_t)nxo * ny * nz);
414 mult(scale);
415#endif //USE_FFTW3
416
417 set_fftpad(true);
418 set_complex(false);
419 if(ny==1 && nz==1) set_complex_x(false);
420 set_ri(false);
421 update();
422
423 EXITFUNC;
424// return this;
425}
426#undef rdata
427
428
429EMBytes EMData::render_ap24(int x0, int y0, int ixsize, int iysize,
430 int bpl, float scale, int mingray, int maxgray,
431 float render_min, float render_max,float gamma,int flags)
432{
433 ENTERFUNC;
434
435 int asrgb;
436 int hist=(flags&2)/2;
437 int invy=(flags&4)?1:0;
438
439 if (!is_complex()) throw ImageDimensionException("complex only");
440
441 if (get_ndim() != 2) {
442 throw ImageDimensionException("2D only");
443 }
444
445 if (is_complex()) ri2ap();
446
447 if (render_max <= render_min) {
448 render_max = render_min + 0.01f;
449 }
450
451 if (gamma<=0) gamma=1.0;
452
453 // Calculating a full floating point gamma for
454 // each pixel in the image slows rendering unacceptably
455 // however, applying a gamma-mapping to an 8 bit colorspace
456 // has unaccepable accuracy. So, we oversample the 8 bit colorspace
457 // as a 12 bit colorspace and apply the gamma mapping to that
458 // This should produce good accuracy for gamma values
459 // larger than 0.5 (and a high upper limit)
460 static int smg0=0,smg1=0; // while this destroys threadsafety in the rendering process
461 static float sgam=0; // it is necessary for speed when rendering large numbers of small images
462 static unsigned char gammamap[4096];
463 if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
464 for (int i=0; i<4096; i++) {
465 if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
466 else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
467 }
468 }
469 smg0=mingray; // so we don't recompute the map unless something changes
470 smg1=maxgray;
471 sgam=gamma;
472
473 if (flags&8) asrgb=4;
474 else if (flags&1) asrgb=3;
475 else throw ImageDimensionException("must set flag 1 or 8");
476
477 EMBytes ret=EMBytes();
478// ret.resize(iysize*bpl);
479 ret.assign(iysize*bpl+hist*1024,char(mingray));
480 unsigned char *data=(unsigned char *)ret.data();
481 unsigned int *histd=(unsigned int *)(data+iysize*bpl);
482 if (hist) {
483 for (int i=0; i<256; i++) histd[i]=0;
484 }
485
486 float rm = render_min;
487 float inv_scale = 1.0f / scale;
488 int ysize = iysize;
489 int xsize = ixsize;
490
491 int ymin = 0;
492 if (iysize * inv_scale > ny) {
493 ymin = (int) (iysize - ny / inv_scale);
494 }
495
496 float gs = (maxgray - mingray) / (render_max - render_min);
497 float gs2 = 4095.999f / (render_max - render_min);
498// float gs2 = 1.0 / (render_max - render_min);
499 if (render_max < render_min) {
500 gs = 0;
501 rm = FLT_MAX;
502 }
503
504 int dsx = -1;
505 int dsy = 0;
506 int remx = 0;
507 int remy = 0;
508 const int scale_n = 100000;
509
510 int addi = 0;
511 int addr = 0;
512 if (inv_scale == floor(inv_scale)) {
513 dsx = (int) inv_scale;
514 dsy = (int) (inv_scale * nx);
515 }
516 else {
517 addi = (int) floor(inv_scale);
518 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
519 }
520
521 int xmin = 0;
522 if (x0 < 0) {
523 xmin = (int) (-x0 / inv_scale);
524 xsize -= (int) floor(x0 / inv_scale);
525 x0 = 0;
526 }
527
528 if ((xsize - xmin) * inv_scale > (nx - x0)) {
529 xsize = (int) ((nx - x0) / inv_scale + xmin);
530 }
531 int ymax = ysize - 1;
532 if (y0 < 0) {
533 ymax = (int) (ysize + y0 / inv_scale - 1);
534 ymin += (int) floor(y0 / inv_scale);
535 y0 = 0;
536 }
537
538 if (xmin < 0) xmin = 0;
539 if (ymin < 0) ymin = 0;
540 if (xsize > ixsize) xsize = ixsize;
541 if (ymax > iysize) ymax = iysize;
542
543 int lmax = nx * ny - 1;
544
545 int mid=nx*ny/2;
546 float* image_data = get_data();
547 if (dsx != -1) {
548 int l = y0 * nx;
549 for (int j = ymax; j >= ymin; j--) {
550 int ll = x0;
551 for (int i = xmin; i < xsize; i++) {
552 if (l + ll > lmax || ll >= nx - 2) break;
553
554 int k = 0;
555 unsigned char p;
556 int ph;
557 if (ll >= nx / 2) {
558 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
559 else k = 2 * (ll - nx / 2) + l + 2 + nx;
560 if (k>=mid) k-=mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
561 else k+=mid;
562 ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
563 }
564 else {
565 k = nx * ny - (l + 2 * ll) - 2;
566 ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
567 if (k>=mid) k-=mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
568 else k+=mid;
569 }
570 float t = image_data[k];
571 if (t <= rm) p = mingray;
572 else if (t >= render_max) p = maxgray;
573 else if (gamma!=1.0) {
574 k=(int)(gs2 * (t-render_min)); // map float value to 0-4096 range
575 p = gammamap[k]; // apply gamma using precomputed gamma map
576 }
577 else {
578 p = (unsigned char) (gs * (t - render_min));
579 p += mingray;
580 }
581 if (ph<256) {
582 data[i * asrgb + j * bpl] = p*(255-ph)/256;
583 data[i * asrgb + j * bpl+1] = p*ph/256;
584 data[i * asrgb + j * bpl+2] = 0;
585 }
586 else if (ph<512) {
587 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
588 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
589 data[i * asrgb + j * bpl] = 0;
590 }
591 else {
592 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
593 data[i * asrgb + j * bpl] = p*(ph-512)/256;
594 data[i * asrgb + j * bpl+1] = 0;
595 }
596 if (hist) histd[p]++;
597 ll += dsx;
598 }
599 l += dsy;
600 }
601 }
602 else {
603 remy = 10;
604 int l = y0 * nx;
605 for (int j = ymax; j >= ymin; j--) {
606 int br = l;
607 remx = 10;
608 int ll = x0;
609 for (int i = xmin; i < xsize - 1; i++) {
610 if (l + ll > lmax || ll >= nx - 2) {
611 break;
612 }
613 int k = 0;
614 unsigned char p;
615 int ph;
616 if (ll >= nx / 2) {
617 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
618 else k = 2 * (ll - nx / 2) + l + 2 + nx;
619 if (k>=mid) k-=mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
620 else k+=mid;
621 ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
622 }
623 else {
624 k = nx * ny - (l + 2 * ll) - 2;
625 if (k>=mid) k-=mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
626 else k+=mid;
627 ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
628 }
629
630 float t = image_data[k];
631 if (t <= rm)
632 p = mingray;
633 else if (t >= render_max) {
634 p = maxgray;
635 }
636 else if (gamma!=1.0) {
637 k=(int)(gs2 * (t-render_min)); // map float value to 0-4096 range
638 p = gammamap[k]; // apply gamma using precomputed gamma map
639 }
640 else {
641 p = (unsigned char) (gs * (t - render_min));
642 p += mingray;
643 }
644 if (ph<256) {
645 data[i * asrgb + j * bpl] = p*(255-ph)/256;
646 data[i * asrgb + j * bpl+1] = p*ph/256;
647 data[i * asrgb + j * bpl+2] = 0;
648 }
649 else if (ph<512) {
650 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
651 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
652 data[i * asrgb + j * bpl] = 0;
653 }
654 else {
655 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
656 data[i * asrgb + j * bpl] = p*(ph-512)/256;
657 data[i * asrgb + j * bpl+1] = 0;
658 }
659 if (hist) histd[p]++;
660 ll += addi;
661 remx += addr;
662 if (remx > scale_n) {
663 remx -= scale_n;
664 ll++;
665 }
666 }
667 l = br + addi * nx;
668 remy += addr;
669 if (remy > scale_n) {
670 remy -= scale_n;
671 l += nx;
672 }
673 }
674 }
675
676 // this replicates r -> g,b
677 if (asrgb==4) {
678 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
679 for (int i=xmin; i<xsize*4; i+=4) {
680 data[i+j+3]=255;
681 }
682 }
683 }
684
685 EXITFUNC;
686
687 // ok, ok, not the most efficient place to do this, but it works
688 if (invy) {
689 int x,y;
690 char swp;
691 for (y=0; y<iysize/2; y++) {
692 for (x=0; x<ixsize; x++) {
693 swp=ret[y*bpl+x];
694 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
695 ret[(iysize-y-1)*bpl+x]=swp;
696 }
697 }
698 }
699
700 // return PyString_FromStringAndSize((const char*) data,iysize*bpl);
701 return ret;
702}
703
704
705void EMData::render_amp24( int x0, int y0, int ixsize, int iysize,
706 int bpl, float scale, int mingray, int maxgray,
707 float render_min, float render_max, void *ref,
708 void cmap(void *, int coord, unsigned char *tri))
709{
710 ENTERFUNC;
711
712 if (get_ndim() != 2) {
713 throw ImageDimensionException("2D only");
714 }
715
716 if (is_complex()) {
717 ri2ap();
718 }
719
720 if (render_max <= render_min) {
721 render_max = render_min + 0.01f;
722 }
723
724 std::string ret=std::string();
725 ret.resize(iysize*bpl);
726 unsigned char *data=(unsigned char *)ret.data();
727
728 float rm = render_min;
729 float inv_scale = 1.0f / scale;
730 int ysize = iysize;
731 int xsize = ixsize;
732 const int scale_n = 100000;
733
734 int ymin = 0;
735 if ( iysize * inv_scale > ny) {
736 ymin = (int) (iysize - ny / inv_scale);
737 }
738 float gs = (maxgray - mingray) / (render_max - render_min);
739 if (render_max < render_min) {
740 gs = 0;
741 rm = FLT_MAX;
742 }
743 int dsx = -1;
744 int dsy = 0;
745 if (inv_scale == floor(inv_scale)) {
746 dsx = (int) inv_scale;
747 dsy = (int) (inv_scale * nx);
748 }
749 int addi = 0;
750 int addr = 0;
751
752 if (dsx == -1) {
753 addi = (int) floor(inv_scale);
754 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
755 }
756
757 int remx = 0;
758 int remy = 0;
759 int xmin = 0;
760 if (x0 < 0) {
761 xmin = (int) (-x0 / inv_scale);
762 xsize -= (int) floor(x0 / inv_scale);
763 x0 = 0;
764 }
765
766 if ((xsize - xmin) * inv_scale > (nx - x0)) {
767 xsize = (int) ((nx - x0) / inv_scale + xmin);
768 }
769 int ymax = ysize - 1;
770 if (y0 < 0) {
771 ymax = (int) (ysize + y0 / inv_scale - 1);
772 ymin += (int) floor(y0 / inv_scale);
773 y0 = 0;
774 }
775
776
777 if (xmin < 0) {
778 xmin = 0;
779 }
780
781 if (ymin < 0) {
782 ymin = 0;
783 }
784 if (xsize > ixsize) {
785 xsize = ixsize;
786 }
787 if (ymax > iysize) {
788 ymax = iysize;
789 }
790
791 int lmax = nx * ny - 1;
792 unsigned char tri[3];
793 float* image_data = get_data();
794 if (is_complex()) {
795 if (dsx != -1) {
796 int l = y0 * nx;
797 for (int j = ymax; j >= ymin; j--) {
798 int ll = x0;
799 for (int i = xmin; i < xsize; i++, ll += dsx) {
800 if (l + ll > lmax || ll >= nx - 2) {
801 break;
802 }
803 int kk = 0;
804 if (ll >= nx / 2) {
805 if (l >= (ny - inv_scale) * nx) {
806 kk = 2 * (ll - nx / 2) + 2;
807 }
808 else {
809 kk = 2 * (ll - nx / 2) + l + 2 + nx;
810 }
811 }
812 else {
813 kk = nx * ny - (l + 2 * ll) - 2;
814 }
815 int k = 0;
816 float t = image_data[kk];
817 if (t <= rm) {
818 k = mingray;
819 }
820 else if (t >= render_max) {
821 k = maxgray;
822 }
823 else {
824 k = (int) (gs * (t - render_min));
825 k += mingray;
826 }
827 tri[0] = static_cast < unsigned char >(k);
828 cmap(ref, kk, tri);
829 data[i * 3 + j * bpl] = tri[0];
830 data[i * 3 + 1 + j * bpl] = tri[1];
831 data[i * 3 + 2 + j * bpl] = tri[2];
832 }
833 l += dsy;
834 }
835 }
836 else {
837 remy = 10;
838 for (int j = ymax, l = y0 * nx; j >= ymin; j--) {
839 int br = l;
840 remx = 10;
841 for (int i = xmin, ll = x0; i < xsize - 1; i++) {
842 if (l + ll > lmax || ll >= nx - 2) {
843 break;
844 }
845 int kk = 0;
846 if (ll >= nx / 2) {
847 if (l >= (ny * nx - nx)) {
848 kk = 2 * (ll - nx / 2) + 2;
849 }
850 else {
851 kk = 2 * (ll - nx / 2) + l + 2 + nx;
852 }
853 }
854 else {
855 kk = nx * ny - (l + 2 * ll) - 2;
856 }
857 int k = 0;
858 float t = image_data[kk];
859 if (t <= rm) {
860 k = mingray;
861 }
862 else if (t >= render_max) {
863 k = maxgray;
864 }
865 else {
866 k = (int) (gs * (t - render_min));
867 k += mingray;
868 }
869 tri[0] = static_cast < unsigned char >(k);
870 cmap(ref, kk, tri);
871 data[i * 3 + j * bpl] = tri[0];
872 data[i * 3 + 1 + j * bpl] = tri[1];
873 data[i * 3 + 2 + j * bpl] = tri[2];
874 ll += addi;
875 remx += addr;
876 if (remx > scale_n) {
877 remx -= scale_n;
878 ll++;
879 }
880 }
881 l = br + addi * nx;
882 remy += addr;
883 if (remy > scale_n) {
884 remy -= scale_n;
885 l += nx;
886 }
887 }
888 }
889 }
890 else {
891 if (dsx != -1) {
892 for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
893 int br = l;
894 for (int i = xmin; i < xsize; i++, l += dsx) {
895 if (l > lmax) {
896 break;
897 }
898 float t = image_data[l];
899 int k = 0;
900 if (t <= rm) {
901 k = mingray;
902 }
903 else if (t >= render_max) {
904 k = maxgray;
905 }
906 else {
907 k = (int) (gs * (t - render_min));
908 k += mingray;
909 }
910 tri[0] = static_cast < unsigned char >(k);
911 cmap(ref, l, tri);
912 data[i * 3 + j * bpl] = tri[0];
913 data[i * 3 + 1 + j * bpl] = tri[1];
914 data[i * 3 + 2 + j * bpl] = tri[2];
915 }
916 l = br + dsy;
917 }
918 }
919 else {
920 remy = 10;
921 for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
922 int br = l;
923 remx = 10;
924 for (int i = xmin; i < xsize; i++) {
925 if (l > lmax) {
926 break;
927 }
928 float t = image_data[l];
929 int k = 0;
930 if (t <= rm) {
931 k = mingray;
932 }
933 else if (t >= render_max) {
934 k = maxgray;
935 }
936 else {
937 k = (int) (gs * (t - render_min));
938 k += mingray;
939 }
940 tri[0] = static_cast < unsigned char >(k);
941 cmap(ref, l, tri);
942 data[i * 3 + j * bpl] = tri[0];
943 data[i * 3 + 1 + j * bpl] = tri[1];
944 data[i * 3 + 2 + j * bpl] = tri[2];
945 l += addi;
946 remx += addr;
947 if (remx > scale_n) {
948 remx -= scale_n;
949 l++;
950 }
951 }
952 l = br + addi * nx;
953 remy += addr;
954 if (remy > scale_n) {
955 remy -= scale_n;
956 l += nx;
957 }
958 }
959 }
960 }
961
962 EXITFUNC;
963}
964
965void EMData::ap2ri()
966{
967 ENTERFUNC;
968
969 if (!is_complex() || is_ri()) {
970 return;
971 }
972
973 Util::ap2ri(get_data(), (size_t)nx * ny * nz);
974 set_ri(true);
975 update();
976 EXITFUNC;
977}
978
979void EMData::ri2inten()
980{
981 ENTERFUNC;
982
983 if (!is_complex()) return;
984 if (!is_ri()) ap2ri();
985
986 float * data = get_data();
987 size_t size = (size_t)nx * ny * nz;
988 for (size_t i = 0; i < size; i += 2) {
989 data[i]=data[i]*data[i]+data[i+1]*data[i+1];
990 data[i+1]=0;
991 }
992
993 set_attr("is_intensity", int(1));
994 update();
995 EXITFUNC;
996}
997
998
999void EMData::ri2ap()
1000{
1001 ENTERFUNC;
1002
1003 if (!is_complex() || !is_ri()) {
1004 return;
1005 }
1006
1007 float * data = get_data();
1008
1009 size_t size = (size_t)nx * ny * nz;
1010 for (size_t i = 0; i < size; i += 2) {
1011 float f = (float)hypot(data[i], data[i + 1]);
1012 if (data[i] == 0 && data[i + 1] == 0) {
1013 data[i + 1] = 0;
1014 }
1015 else {
1016 data[i + 1] = atan2(data[i + 1], data[i]);
1017 }
1018 data[i] = f;
1019 }
1020
1021 set_ri(false);
1022 update();
1023 EXITFUNC;
1024}
1025
1026
1027float calc_bessel(const int n, const float& x) {
1028 gsl_sf_result result;
1029// int success =
1030 gsl_sf_bessel_Jn_e(n,(double)x, &result);
1031 return (float)result.val;
1032}
1033
1034EMData* EMData::bispecRotTransInvN(int N, int NK)
1035{
1036
1037 int EndP = this -> get_xsize(); // length(fTrueVec);
1038 int Mid = (int) ((1+EndP)/2);
1039 int End = 2*Mid-1;
1040
1041 int CountxyMax = End*End;
1042
1043 int *SortfkInds = new int[CountxyMax];
1044 int *kVecX = new int[CountxyMax];
1045 int *kVecY = new int[CountxyMax];
1046 float *fkVecR = new float[CountxyMax];
1047 float *fkVecI = new float[CountxyMax];
1048 float *absD1fkVec = new float[CountxyMax];
1049 float *absD1fkVecSorted = new float[CountxyMax];
1050
1051 float *jxjyatan2 = new float[End*End]; // jxjyatan2[jy*End + jx] = atan2(jy+1-Mid , jx +1 -Mid);
1052
1053 EMData * ThisCopy = new EMData(End,End);
1054
1055 for (int jx=0; jx <End ; jx++) {
1056 for (int jy=0; jy <End ; jy++) {
1057 float ValNow = this -> get_value_at(jx,jy);
1058 ThisCopy -> set_value_at(jx,jy,ValNow);
1059// cout<< " jxM= " << jx+1<<" jyM= " << jy+1<< "ValNow" << ValNow << endl; // Works
1060 }}
1061
1062
1063 EMData* fk = ThisCopy -> do_fft();
1064 fk ->process_inplace("xform.fourierorigin.tocenter");
1065
1066// EMData* fk
1067 EMData* fkRCopy = new EMData(End,End);
1068 EMData* fkICopy = new EMData(End,End);
1069 EMData* fkCopy = new EMData(End,End);
1070
1071
1072 for (int jCount= 0; jCount<End*End; jCount++) {
1073// jCount = jy*End + jx;
1074 int jx = jCount%End ;
1075 int jy = (jCount-jx)/End ;
1076 jxjyatan2[jCount] = atan2((float)(jy+1-Mid) , (float)(jx +1-Mid));
1077 }
1078
1079
1080 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) { // kEx twice the value of the Fourier
1081 // x variable: EMAN index for real, imag
1082 int kx = kEx/2; // kx is the value of the Fourier variable
1083 int kIx = kx+Mid-1; // This is the value of the index for a matlab image (-1)
1084 int kCx = -kx ;
1085 int kCIx = kCx+ Mid-1 ;
1086 for (int kEy= 0 ; kEy<End; kEy++) { // This is the value of the EMAN index
1087 int kIy = kEy ; // This is the value of the index for a matlab image (-1)
1088 int ky = kEy+1-Mid; // (kEy+ Mid-1)%End - Mid+1 ; // This is the actual value of the Fourier variable
1089 float realVal = fk -> get_value_at(kEx ,kEy) ;
1090 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
1091 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
1092 float fkAng = atan2(imagVal,realVal);
1093
1094 float NewRealVal ;
1095 float NewImagVal ;
1096 float AngMatlab ;
1097
1098 if (kIx==Mid-1) {
1099// AngMatlab = -fkAng - 2.*M_PI*(kIy+ 1-Mid)*(Mid)/End;
1100// cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]=" << fkVecI[i] <<" angle[i]= " << AngMatlab << endl;
1101 }
1102
1103 if (kIx>Mid-1){
1104// cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]=" << fkVecI[i] <<" angle[i]= " << AngMatlab << endl;
1105 }
1106
1107 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
1108 NewRealVal = absVal*cos(AngMatlab);
1109 NewImagVal = absVal*sin(AngMatlab);
1110
1111
1112 fkVecR[kIy+kIx *End] = NewRealVal ;
1113 fkVecR[kIy+kCIx*End] = NewRealVal ;
1114 fkVecI[kIy+kIx *End] = NewImagVal ;
1115 fkVecI[kIy+kCIx*End] = -NewImagVal ;
1116 absD1fkVec[kIy + kIx *End] = absVal;
1117 absD1fkVec[kIy + kCIx *End] = absVal;
1118 kVecX[kIy+kIx *End] = kx ;
1119 kVecX[kIy+kCIx *End] = kCx ;
1120 kVecY[kIy+kIx *End] = ky ;
1121 kVecY[kIy+kCIx *End] = ky ;
1122// printf("kx=%d,ky=%d,tempVal =%f+ i %4.2f \n",kx,ky,realVal,imagVal );
1123// cout << "kx = " << kx << "; ky = "<< ky << "; val is" << realVal<<"+ i "<<imagVal<< endl;
1124
1125// cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; fkAng*9/ 2pi is " << fkAng*9/2/M_PI<< endl;
1126// cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; absval is " << absVal<< "; realval is " << NewRealVal<< "; imagval is " << NewImagVal<< endl;
1127 fkCopy -> set_value_at(kIx ,kIy, absVal);
1128 fkCopy -> set_value_at(kCIx,kIy, absVal);
1129 fkRCopy -> set_value_at(kIx, kIy, NewRealVal);
1130 fkRCopy -> set_value_at(kCIx,kIy, NewRealVal);
1131 fkICopy -> set_value_at(kIx, kIy, NewImagVal);
1132 fkICopy -> set_value_at(kCIx,kIy,-NewImagVal);
1133
1134 }
1135 }
1136#ifdef _WIN32
1137 _unlink("fkCopy.???");
1138 _unlink("fk?Copy.???");
1139#else
1140 int rslt;
1141 rslt = system("rm -f fkCopy.???"); rslt++;
1142 rslt = system("rm -f fk?Copy.???"); rslt++;
1143#endif //_WIN32
1144 fkCopy -> write_image("fkCopy.img");
1145 fkRCopy -> write_image("fkRCopy.img");
1146 fkICopy -> write_image("fkICopy.img");
1147
1148 cout << "Starting the sort "<< endl;
1149
1150 vector< pair<float, int> > absInds;
1151 for(int i = 0; i < CountxyMax; ++i ) {
1152 pair<float,int> p;
1153 p = make_pair(absD1fkVec[i],i); // p = make_pair(rand(),i);
1154 absInds.push_back( p);
1155 }
1156
1157 std::sort(absInds.begin(),absInds.end());
1158
1159 for(int i = 0; i < CountxyMax; ++i ) {
1160 pair<float,int> p ;
1161 p = absInds[i] ;
1162 absD1fkVecSorted[CountxyMax-1-i] = p.first ;
1163 SortfkInds[CountxyMax-1-i] = p.second ;
1164 }
1165
1166 cout << "Ending the sort "<< endl;
1167
1168// 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};
1169
1170
1171 for(int i = 0; i < CountxyMax; ++i ) { // creates a new fkVec
1172// int Si = SortfkInds[i];
1173// int kIx = (int) Si/End; kIx = (int) i/End; // i = kIx*End+kIy
1174// int kIy = Si - kIx*End; kIy = i - kIx*End;
1175// int iC = (End-1-kIx)*End + (End-1-kIy);
1176// if (i<30) { cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " valAft=" << absD1fkVecSorted[i]<< " valBef=" << absD1fkVec[Si] << " SortfkInds = " << Si <<endl; }// This worked
1177// cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]=" << fkVecI[i] <<" angle[i]= " << fkAng << endl;
1178 }
1179 cout<< "Ratio of Last Amplitude to First Amplitude= " << absD1fkVecSorted[NK] /absD1fkVecSorted[0] << endl;
1180
1181// pause;
1182
1183// for(int i = 0; i < NK; ++i ) { // Prints out the new fkVec , CountxyMax
1184// int Si= SortfkInds[i];
1185// int kIx = (int) Si/End; // i = kIx*End+kIy
1186// int kIy = Si - kIx*End;
1187// 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;
1188// }
1189
1190// angEMAN+angMat+angDiff =0 mod 2 pi
1191
1192// angDiff= 2*pi*((-4):4)*(Mid)/End; angEMAN+angMat+angDiff= integer*2 *pi
1193// [ absD1fkVecSorted, SortfkInds] =sort( absD1fkVec,'descend') ;
1194// Util::sort_mat(&absD1fkVec[0],&absD1fkVec[Countxy],&SortfkInds[0],&SortfkInds[Countxy]);
1195
1196
1197// Let radial sampling be 0:0.5:(Mid-1)
1198
1199 // int NK= min(12,CountxyMax) ;
1200
1201
1202
1203 cout << "NK = " << NK << endl;
1204 float frR= 3.0/4.0;
1205 int LradRange= (int) (floor(Mid/frR)) ;
1206
1207 float *radRange = new float[LradRange]; //= 0:.75:(Mid-1);
1208 radRange[0]=0;
1209 for (int irad=1; irad < LradRange; irad++){
1210 radRange[irad] = radRange[irad-1] + frR; }
1211
1212
1213
1214 // should equal to (2*Mid-1)
1215 cout << "Starting the calculation of invariants for N= " << N << endl;
1216
1217/* int NMax=5; */
1218
1219 EMData* RotTransInv = new EMData();
1220 RotTransInv -> set_size(LradRange,LradRange);
1221
1222
1223// float *RotTransInv = new float[LradRange*LradRange ] ;
1224// float *RotTransInvN = new float[LradRange*LradRange*(NMax+1) ] ;
1225
1226// for (int N=0 ; N<NMax; N++) {
1227
1228 for (int jr1=0; jr1 < LradRange ; jr1++ ) { // LradRange
1229 float r1= radRange[jr1];
1230// cout << "Pre jr2 "<< endl;
1231 for (int jr2=0; jr2<LradRange; jr2++ ) { //LradRange
1232 float r2= radRange[jr2];
1233 float RotTransInvTemp=0;
1234 for (int jCountkxy =0; jCountkxy<NK; jCountkxy++){
1235 int Countkxy =SortfkInds[jCountkxy] ; // 1: CountxyMax
1236 int kx = kVecX[Countkxy] ;
1237 int ky = kVecY[Countkxy] ;
1238 float k2 = (float)(kx*kx+ky*ky);
1239 if (k2==0) { continue;}
1240 float phiK =0;
1241 if (k2>0) phiK= jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1]; phiK= atan2((float)ky,(float)kx);
1242
1243 float fkR = fkVecR[Countkxy] ;
1244 float fkI = fkVecI[Countkxy] ;
1245/* printf("jCountkxy=%d, Countkxy=%d,absD1fkVec(Countkxy)=%f,\t\t kx=%d, ky=%d \n", jCountkxy, Countkxy, absD1fkVec[Countkxy], kx, ky);*/
1246
1247 for (int jCountqxy =0; jCountqxy<NK; jCountqxy++){
1248 int Countqxy =SortfkInds[jCountqxy] ; // Countqxy is the index for absD1fkVec
1249 int qx = kVecX[Countqxy] ;
1250 int qy = kVecY[Countqxy] ;
1251 int q2 = qx*qx+qy*qy;
1252 if (q2==0) {continue;}
1253 float phiQ =0;
1254 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1]; phiQ=atan2((float)qy,(float)qx);
1255 float fqR = fkVecR[Countqxy] ;
1256 float fqI = fkVecI[Countqxy] ;
1257 int kCx = (-kx-qx);
1258 int kCy = (-ky-qy);
1259 int kCIx = ((kCx+Mid+2*End)%End);// labels of the image in C
1260 int kCIy = ((kCy+Mid+2*End)%End);
1261 kCx = kCIx-Mid; // correct
1262 kCy = kCIy-Mid; // correct
1263 int CountCxy = kCIx*End+kCIy;
1264 float fCR = fkVecR[CountCxy];
1265 float fCI = fkVecI[CountCxy];
1266 if (jr1+jr2==-1) {
1267 printf("jCountqxy=%d , Countqxy=%d, absD1fkVec(Countqxy)=%f,qx=%d, qy=%d \n", jCountqxy, Countqxy, absD1fkVec[Countqxy],qx, qy);
1268 printf(" CountCxy=%d,absD1fkVec[CountCxy]=%f, kCx=%d, kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
1269 }
1270 for (int p=0; p<NK; p++){
1271// printf("p=%d, SortfkInds[p]=%d, CountCxy =%d \n", p,SortfkInds[p], CountCxy);
1272 if (SortfkInds[p]==CountCxy){
1273 float Arg1 = 2.0f*M_PI*r1*::sqrt((float) q2)/End;
1274 float Arg2 = 2.0f*M_PI*r2*::sqrt((float) k2)/End;
1275// printf("Arg1=%4.2f, Arg2=%4.2f, \n",Arg1, Arg2 );
1276// if (Arg1+ Arg2<15) {
1277 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI))
1278 * cos(N*(phiK-phiQ+M_PI));
1279 bispectemp -= (fkR*(fqR*fCI + fqI*fCR) +fkI*(fqR*fCR - fqI*fCI))
1280 * sin(N*(phiK-phiQ+M_PI));
1281 float bess1 = calc_bessel(N, Arg1 );
1282 float bess2 = calc_bessel(N, Arg2 );
1283// printf("fkr=%4.2f, fqr=%4.2f, bess1=%4.2f,bess2=%4.2f \n",fkR, fqR, bess1, bess2);
1284/* printf("p =%d, SortfkInds[p]=%d, CountCxy=%d, Arg1 =%4.2f, bess1=%4.2f, \n",
1285 p, SortfkInds[p],CountCxy, Arg1, bess1);*/
1286 RotTransInvTemp = RotTransInvTemp + bispectemp * bess1*bess2 ;
1287// }
1288 }
1289 }
1290 } // jCountqxy
1291 } // jCountkxy
1292 RotTransInv -> set_value_at(jr1,jr2, RotTransInvTemp) ;
1293/* RotTransInvN[jr1 + LradRange*jr2+LradRange*LradRange*N] = RotTransInvTemp ;*/
1294 } //jr2
1295 } //jr1
1296// }//N
1297
1298 return RotTransInv ;
1299
1300
1301}
1302
1303
1304
1305/*
1306// find example
1307#include <iostream>
1308#include <algorithm>
1309#include <vector>
1310using namespace std;
1311
1312int main () {
1313 int myints[] = { 10, 20, 30 ,40 };
1314 int * p;
1315
1316 // pointer to array element:
1317 p = find(myints,myints+4,30);
1318 ++p;
1319 cout << "The element following 30 is " << *p << endl;
1320
1321 vector<int> myvector (myints,myints+4);
1322 vector<int>::iterator it;
1323
1324 // iterator to vector element:
1325 it = find (myvector.begin(), myvector.end(), 30);
1326 ++it;
1327 cout << "The element following 30 is " << *it << endl;
1328
1329 return 0;
1330}*/
1331
1332EMData* EMData::bispecRotTransInvDirect(int type)
1333{
1334
1335 int EndP = this -> get_xsize(); // length(fTrueVec);
1336 int Mid = (int) ((1+EndP)/2);
1337 int End = 2*Mid-1;
1338
1339 int CountxyMax = End*End;
1340
1341// int *SortfkInds = new int[CountxyMax];
1342 int *kVecX = new int[CountxyMax];
1343 int *kVecY = new int[CountxyMax];
1344 float *fkVecR = new float[CountxyMax];
1345 float *fkVecI = new float[CountxyMax];
1346 float *absD1fkVec = new float[CountxyMax];
1347// float *absD1fkVecSorted = new float[CountxyMax];
1348
1349
1350 float *jxjyatan2 = new float[End*End];
1351
1352
1353 EMData * ThisCopy = new EMData(End,End);
1354
1355 for (int jx=0; jx <End ; jx++) { // create jxjyatan2
1356 for (int jy=0; jy <End ; jy++) {
1357 float ValNow = this -> get_value_at(jx,jy);
1358 ThisCopy -> set_value_at(jx,jy,ValNow);
1359 jxjyatan2[jy*End + jx] = atan2((float)(jy+1-Mid) , (float)(jx +1 -Mid));
1360// cout<< " jxM= " << jx+1<<" jyM= " << jy+1<< "ValNow" << ValNow << endl; // Works
1361 }}
1362
1363
1364 EMData* fk = ThisCopy -> do_fft();
1365 fk ->process_inplace("xform.fourierorigin.tocenter");
1366
1367// Create kVecX , kVecy etc
1368
1369 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) { // kEx twice the value of the Fourier
1370 // x variable: EMAN index for real, imag
1371 int kx = kEx/2; // kx is the value of the Fourier variable
1372 int kIx = kx+Mid-1; // This is the value of the index for a matlab image (-1)
1373 int kCx = -kx ;
1374 int kCIx = kCx+ Mid-1 ;
1375 for (int kEy= 0 ; kEy<End; kEy++) { // This is the value of the EMAN index
1376 int kIy = kEy ; // This is the value of the index for a matlab image (-1)
1377 int ky = kEy+1-Mid; // (kEy+ Mid-1)%End - Mid+1 ; // This is the actual value of the Fourier variable
1378 float realVal = fk -> get_value_at(kEx ,kEy) ;
1379 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
1380 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
1381 float fkAng = atan2(imagVal,realVal);
1382
1383 float NewRealVal ;
1384 float NewImagVal ;
1385 float AngMatlab ;
1386
1387 if (kIx==Mid-1) {
1388// AngMatlab = -fkAng - 2.*M_PI*(kIy+ 1-Mid)*(Mid)/End;
1389 }
1390
1391 if (kIx>Mid-1){
1392// cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]=" << fkVecI[i] <<" angle[i]= " << AngMatlab << endl;
1393 }
1394
1395 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
1396 NewRealVal = absVal*cos(AngMatlab);
1397 NewImagVal = absVal*sin(AngMatlab);
1398
1399
1400 fkVecR[ kIy +kIx *End] = NewRealVal ;
1401 fkVecR[(End-1-kIy)+kCIx*End] = NewRealVal ;
1402 fkVecI[ kIy +kIx *End] = NewImagVal ;
1403 fkVecI[(End-1-kIy)+kCIx*End] = -NewImagVal ;
1404 absD1fkVec[(End-1-kIy) + kIx *End] = absVal;
1405 absD1fkVec[(End-1-kIy) + kCIx *End] = absVal;
1406 kVecX[kIy+kIx *End] = kx ;
1407 kVecX[kIy+kCIx *End] = kCx ;
1408 kVecY[kIy+kIx *End] = ky ;
1409 kVecY[kIy+kCIx *End] = ky ;
1410
1411 // cout << " kIxM= " << kIx+1 << " kIy=" << kIy+1 << " fkVecR[i] =" << NewRealVal << " fkVecI[i]=" << NewImagVal <<" angle[i]= " << AngMatlab << " Total Index" << kIy+kIx *End << endl;
1412
1413// printf("kx=%d,ky=%d,tempVal =%f+ i %4.2f \n",kx,ky,realVal,imagVal );
1414// cout << "kx = " << kx << "; ky = "<< ky << "; val is" << realVal<<"+ i "<<imagVal<< endl;
1415
1416// cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; fkAng*9/ 2pi is " << fkAng*9/2/M_PI<< endl;
1417// cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; absval is " << absVal<< "; realval is " << NewRealVal<< "; imagval is " << NewImagVal<< endl;
1418 }
1419 }
1420
1421
1422// for (int TotalInd = 0 ; TotalInd < CountxyMax ; TotalInd++){
1423// int kx = kVecX[TotalInd]; // This is the value of the index for a matlab image (-1)
1424// int kIx = kx+Mid-1; // This is the value of the index for a matlab image (-1)
1425// int ky = kVecY[TotalInd];
1426// int kIy = ky+Mid-1; // This is the value of the index for a matlab image (-1)
1427 //float fkR = fkVecR[kIy+kIx *End] ;
1428 //float fkI = fkVecI[kIy+kIx *End] ;
1429// }
1430
1431 float frR= 3.0/4.0;
1432 frR= 1;
1433 int LradRange= (int) (1+floor(Mid/frR -.1)) ;
1434
1435 float *radRange = new float[LradRange]; //= 0:.75:(Mid-1);
1436 for (int irad=0; irad < LradRange; irad++){
1437 radRange[irad] = frR*irad;
1438// cout << " irad = " << irad << " radRange[irad]= " << radRange[irad] << " LradRange= " << LradRange << endl;
1439 }
1440
1441 cout << "Starting the calculation of invariants" << endl;
1442
1443
1444 if (type==0) {
1445 int LthetaRange = 59;
1446 float ftR = (2.0f*M_PI/LthetaRange );
1447 float *thetaRange = new float[LthetaRange]; //= 0:.75:(Mid-1);
1448
1449 for (int ith=0; ith < LthetaRange; ith++){
1450 thetaRange[ith] = ftR*ith; }
1451
1452 int TotalVol = LradRange*LradRange*LthetaRange;
1453
1454 float *RotTransInv = new float[TotalVol];
1455 float *WeightInv = new float[TotalVol];
1456
1457 for (int jW=0; jW<TotalVol; jW++) {
1458 RotTransInv[jW] = 0;
1459 WeightInv[jW] = 0;
1460 }
1461
1462 for (int jW=0; jW<TotalVol; jW++) {
1463 RotTransInv[jW] = 0;
1464 WeightInv[jW] = 0;
1465 }
1466 // float *RotTransInv = new float[LradRange*LradRange ] ;
1467 // float *RotTransInvN = new float[LradRange*LradRange*(NMax+1) ] ;
1468
1469 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){ // Main Section for type 0
1470 int kx = kVecX[Countkxy] ;
1471 int ky = kVecY[Countkxy] ;
1472 float k2 = ::sqrt((float)(kx*kx+ky*ky));
1473 float phiK =0;
1474 if (k2>0) phiK = jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1]; // phiK=atan2(ky,kx);
1475 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
1476 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
1477 // printf("Countkxy=%d,\t kx=%d, ky=%d, fkR=%3.2f,fkI=%3.2f \n", Countkxy, kx, ky, fkR, fkI);
1478
1479 if ((k2==0)|| (k2>Mid) ) { continue;}
1480
1481 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){ // This is the innermost loop
1482 int qx = kVecX[Countqxy] ;
1483 int qy = kVecY[Countqxy] ;
1484 float q2 = ::sqrt((float)(qx*qx+qy*qy));
1485 if ((q2==0)|| (q2>Mid) ) {continue;}
1486 float phiQ =0;
1487 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1]; // phiQ=atan2(qy,qx);
1488 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
1489 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
1490 int kCx = (-kx-qx);
1491 int kCy = (-ky-qy);
1492 int kCIx = ((kCx+Mid+2*End)%End);// labels of the image in C
1493 int kCIy = ((kCy+Mid+2*End)%End);
1494 kCx = ((kCIx+End-1)%End)+1-Mid; // correct
1495 kCy = ((kCIy+End-1)%End)+1-Mid ; // correct
1496
1497// float C2 = ::sqrt((float)(kCx*kCx+ kCy*kCy));
1498 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
1499 float fCR = fkVecR[CountCxy];
1500 float fCI = fkVecI[CountCxy];
1501 /* if (Countkxy==1) {
1502 printf(" Countqxy=%d, absD1fkVec(Countqxy)=%f,qx=%d, qy=%d \n", Countqxy, absD1fkVec[Countqxy],qx, qy);
1503 printf(" CountCxy=%d, absD1fkVec[CountCxy]=%f,kCx=%d,kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
1504 }*/
1505// float phiC = jxjyatan2[ (kCy+Mid-1)*End + kCx+Mid-1];
1506 float phiQK = (4*M_PI+phiQ-phiK);
1507 while (phiQK> (2*M_PI)) phiQK -= (2*M_PI);
1508
1509
1510
1511 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
1512
1513 if ((q2<k2) ) continue;
1514// if ((q2<k2) || (C2<k2) || (C2<q2)) continue;
1515
1516 // printf(" CountCxy=%d, absD1fkVec[CountCxy]=%f,kCx=%d,kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
1517
1518 // up to here, matched perfectly with Matlab
1519
1520 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
1521 int k2IndHi = k2IndLo;
1522 float k2Lo= radRange[k2IndLo];
1523 if (k2IndLo+1< LradRange) {
1524 k2IndHi = k2IndLo+1;
1525 }
1526// float k2Hi= radRange[k2IndHi];
1527
1528 float kCof =k2-k2Lo;
1529
1530 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
1531 int q2IndHi=q2IndLo;
1532 float q2Lo= radRange[q2IndLo];
1533 if (q2IndLo+1 < LradRange) {
1534 q2IndHi = q2IndLo+1 ;
1535 }
1536 float qCof = q2-q2Lo;
1537
1538 if ((qCof<0) || (qCof >1) ) {
1539 cout<< "Weird! qCof="<< qCof << " q2="<< q2 << " q2IndLo="<< q2IndLo << endl ;
1540 int x ;
1541 cin >> x ;
1542 }
1543
1544 int thetaIndLo = 0; while ((phiQK>=thetaRange[thetaIndLo+1])&& (thetaIndLo+1<LthetaRange)) thetaIndLo +=1;
1545 int thetaIndHi = thetaIndLo;
1546
1547 float thetaLo = thetaRange[thetaIndLo];
1548 float thetaHi = thetaLo;
1549 float thetaCof = 0;
1550
1551 if (thetaIndLo+1< LthetaRange) {
1552 thetaIndHi = thetaIndLo +1;
1553 }else{
1554 thetaIndHi=0;
1555 }
1556
1557 thetaHi = thetaRange[thetaIndHi];
1558
1559 if (thetaHi==thetaLo) {
1560 thetaCof =0 ;
1561 } else {
1562 thetaCof = (phiQK-thetaLo)/(thetaHi-thetaLo);
1563 }
1564
1565 if ((thetaCof>2*M_PI) ) {
1566 cout<< "Weird! thetaCof="<< thetaCof <<endl ;
1567 thetaCof=0;
1568 }
1569
1570
1571 // if ((thetaIndLo>=58) || (k2IndLo >= LradRange-1) || (q2IndLo >= LradRange-1) ) {
1572
1573
1574 for (int jk =1; jk<=2; jk++){
1575 for (int jq =1; jq<=2; jq++){
1576 for (int jtheta =1; jtheta<=2; jtheta++){
1577
1578 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1))
1579 * (thetaCof+(1-2*thetaCof)*(jtheta==1));
1580
1581
1582 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
1583 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
1584 int thetaInd = thetaIndLo*(jtheta==1) + thetaIndHi*(jtheta ==2);
1585 int TotalInd = thetaInd*LradRange*LradRange+q2Ind*LradRange+k2Ind;
1586 /* if (TotalInd+1 >= LthetaRange*LradRange*LradRange) {
1587 cout << "Weird!!! TotalInd="<< TotalInd << " IndMax" << LthetaRange*LradRange*LradRange << " LradRange=" << LradRange << endl;
1588 cout << "k2Ind= "<< k2Ind << " q2Ind="<< q2Ind << " thetaInd="<< thetaInd << " q2IndLo="<< q2IndLo << " q2IndHi="<< q2IndHi << endl;
1589 cout << "k2=" << k2 << "q2=" << q2 << " phiQK=" << phiQK*180.0/M_PI<< endl;
1590 }*/
1591
1592 RotTransInv[TotalInd] += Weight*bispectemp;
1593 WeightInv[TotalInd] += Weight;
1594 // cout << "k2Ind= "<< k2Ind << " q2Ind="<< q2Ind << "Weight=" << Weight << endl;
1595 }}}
1596 } // Countqxy
1597 } // Countkxy
1598
1599 cout << "Finished Main Section " << endl;
1600
1601 /* RotTransInvN[jr1 + LradRange*jr2+LradRange*LradRange*N] = RotTransInvTemp ;*/
1602
1603 cout << " LradRange " <<LradRange <<" LthetaRange " << LthetaRange << endl;
1604 EMData *RotTransInvF = new EMData(LradRange,LradRange,LthetaRange);
1605 EMData *WeightImage = new EMData(LradRange,LradRange,LthetaRange);
1606
1607 // cout << "FFFFFFF" << endl;
1608 //
1609 // RotTransInvF -> set_size(LradRange,LradRange,LthetaRange);
1610 //
1611 // cout << "GGGG" << endl;
1612
1613 for (int jtheta =0; jtheta < LthetaRange; jtheta++){ // write out to RotTransInvF
1614 for (int jq =0; jq<LradRange; jq++){ // LradRange
1615 for (int jk =0; jk<LradRange ; jk++){// LradRange
1616 // cout << "Hi There" << endl;
1617 int TotalInd = jtheta*LradRange*LradRange+jq*LradRange+jk;
1618 float Weight = WeightInv[TotalInd];
1619 WeightImage -> set_value_at(jk,jq,jtheta,Weight);
1620 RotTransInvF -> set_value_at(jk,jq,jtheta,0);
1621 if (Weight <= 0) continue;
1622 RotTransInvF -> set_value_at(jk,jq,jtheta,RotTransInv[TotalInd] / Weight);// include /Weight
1623 int newjtheta = (LthetaRange- jtheta)%LthetaRange;
1624 RotTransInvF -> set_value_at(jq,jk,newjtheta,RotTransInv[TotalInd]/Weight );// include /Weight
1625 }
1626 }
1627 }
1628
1629 cout << " Almost Done " << endl;
1630#ifdef _WIN32
1631 _unlink("WeightImage.???");
1632#else
1633 int rslt;
1634 rslt = system("rm -f WeightImage.???"); rslt++;
1635#endif //_WIN32
1636 WeightImage -> write_image("WeightImage.img");
1637
1638 return RotTransInvF ;
1639 } // Finish type 0
1640
1641 if (type==1) {
1642 int TotalVol = LradRange*LradRange;
1643
1644 float *RotTransInv = new float[TotalVol];
1645 float *WeightInv = new float[TotalVol];
1646
1647 for (int jW=0; jW<TotalVol; jW++) {
1648 RotTransInv[jW] = 0;
1649 WeightInv[jW] = 0;
1650 }
1651
1652
1653 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){
1654 int kx = kVecX[Countkxy] ;
1655 int ky = kVecY[Countkxy] ;
1656 float k2 = ::sqrt((float)(kx*kx+ky*ky));
1657 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
1658 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
1659 // printf("Countkxy=%d,\t kx=%d, ky=%d, fkR=%3.2f,fkI=%3.2f \n", Countkxy, kx, ky, fkR, fkI);
1660
1661 if ((k2==0)|| (k2>Mid) ) { continue;}
1662
1663 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){ // This is the innermost loop
1664
1665// up to here, matched perfectly with Matlab
1666 int qx = kVecX[Countqxy] ;
1667 int qy = kVecY[Countqxy] ;
1668 float q2 = ::sqrt((float)(qx*qx+qy*qy));
1669 if ((q2==0)|| (q2>Mid) ) {continue;}
1670 if ((q2<k2) ) continue;
1671
1672 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
1673 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
1674
1675 int kCx = (-kx-qx);
1676 int kCy = (-ky-qy);
1677 int kCIx = ((kCx+Mid+2*End)%End);// labels of the image in C
1678 int kCIy = ((kCy+Mid+2*End)%End);
1679 kCx = ((kCIx+End-1)%End)+1-Mid; // correct
1680 kCy = ((kCIy+End-1)%End)+1-Mid ; // correct
1681
1682// float C2 = ::sqrt((float)(kCx*kCx+ kCy*kCy));
1683 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
1684 float fCR = fkVecR[CountCxy];
1685 float fCI = fkVecI[CountCxy];
1686
1687
1688 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
1689
1690
1691 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
1692 int k2IndHi = k2IndLo;
1693 float k2Lo= radRange[k2IndLo];
1694 if (k2IndLo+1< LradRange) {
1695 k2IndHi = k2IndLo+1;
1696 }
1697// float k2Hi= radRange[k2IndHi];
1698
1699 float kCof =k2-k2Lo;
1700
1701 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
1702 int q2IndHi=q2IndLo;
1703 float q2Lo= radRange[q2IndLo];
1704 if (q2IndLo+1 < LradRange) {
1705 q2IndHi = q2IndLo+1 ;
1706 }
1707 float qCof = q2-q2Lo;
1708
1709
1710 for (int jk =1; jk<=2; jk++){
1711 for (int jq =1; jq<=2; jq++){
1712
1713 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1));
1714
1715 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
1716 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
1717 int TotalInd = q2Ind*LradRange+k2Ind;
1718 RotTransInv[TotalInd] += Weight*bispectemp;
1719 WeightInv[TotalInd] += Weight;
1720 // cout << "k2Ind= "<< k2Ind << " q2Ind="<< q2Ind << "Weight=" << Weight << endl;
1721 }}
1722 } // Countqxy
1723 } // Countkxy
1724
1725// cout << "Finished Main Section " << endl;
1726// cout << " LradRange " <<LradRange << endl;
1727
1728
1729 EMData *RotTransInvF = new EMData(LradRange,LradRange);
1730 EMData *WeightImage = new EMData(LradRange,LradRange);
1731
1732 for (int jk =0; jk<LradRange ; jk++){// LradRange
1733 for (int jq =jk; jq<LradRange; jq++){ // LradRange
1734 int TotalInd = jq*LradRange+jk;
1735 int TotalIndBar = jq*LradRange+jk;
1736 float Weight = WeightInv[TotalInd] + WeightInv[TotalIndBar];
1737 if (Weight <=0) continue;
1738 WeightImage -> set_value_at(jk,jq,Weight);
1739 WeightImage -> set_value_at(jq,jk,Weight);
1740 float ValNow = cbrt( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight ) ;
1741 RotTransInvF -> set_value_at(jk,jq,ValNow);// include /Weight
1742 RotTransInvF -> set_value_at(jq,jk,ValNow );// include /Weight
1743 }}
1744
1745#ifdef _WIN32
1746 _unlink("WeightImage.???");
1747#else
1748 int rslt;
1749 rslt = system("rm -f WeightImage.???"); rslt++;
1750#endif //_WIN32
1751 WeightImage -> write_image("WeightImage.img");
1752
1753 return RotTransInvF ;
1754 }
1755 return 0;
1756}
1757
1758
1759void EMData::insert_clip(const EMData * const block, const IntPoint &origin) {
1760 int nx1 = block->get_xsize();
1761 int ny1 = block->get_ysize();
1762 int nz1 = block->get_zsize();
1763
1764 Region area(origin[0], origin[1], origin[2],nx1, ny1, nz1);
1765
1766 //make sure the block fits in EMData
1767 int x0 = (int) area.origin[0];
1768 x0 = x0 < 0 ? 0 : x0;
1769
1770 int y0 = (int) area.origin[1];
1771 y0 = y0 < 0 ? 0 : y0;
1772
1773 int z0 = (int) area.origin[2];
1774 z0 = z0 < 0 ? 0 : z0;
1775
1776 int x1 = (int) (area.origin[0] + area.size[0]);
1777 x1 = x1 > nx ? nx : x1;
1778
1779 int y1 = (int) (area.origin[1] + area.size[1]);
1780 y1 = y1 > ny ? ny : y1;
1781
1782 int z1 = (int) (area.origin[2] + area.size[2]);
1783 z1 = z1 > nz ? nz : z1;
1784 if (z1 <= 0) {
1785 z1 = 1;
1786 }
1787
1788 int xd0 = (int) (area.origin[0] < 0 ? -area.origin[0] : 0);
1789 int yd0 = (int) (area.origin[1] < 0 ? -area.origin[1] : 0);
1790 int zd0 = (int) (area.origin[2] < 0 ? -area.origin[2] : 0);
1791
1792 if (x1 < x0 || y1 < y0 || z1 < z0) return; // out of bounds, this is fine, nothing happens
1793
1794 size_t clipped_row_size = (x1-x0) * sizeof(float);
1795 size_t src_secsize = (size_t)(nx1 * ny1);
1796 size_t dst_secsize = (size_t)(nx * ny);
1797
1798/*
1799#ifdef EMAN2_USING_CUDA
1800 if(block->cudarwdata){
1801 // this is VERY slow.....
1802 float *cudasrc = block->cudarwdata + zd0 * src_secsize + yd0 * nx1 + xd0;
1803 if(!cudarwdata) rw_alloc();
1804 float *cudadst = cudarwdata + z0 * dst_secsize + y0 * nx + x0;
1805 for (int i = z0; i < z1; i++) {
1806 for (int j = y0; j < y1; j++) {
1807 //printf("%x %x %d\n", cudadst, cudasrc, clipped_row_size);
1808 cudaMemcpy(cudadst,cudasrc,clipped_row_size,cudaMemcpyDeviceToDevice);
1809 cudasrc += nx1;
1810 cudadst += nx;
1811 }
1812 cudasrc += src_gap;
1813 cudadst += dst_gap;
1814 }
1815 return;
1816 }
1817#endif
1818*/
1819 float *src = block->get_data() + (size_t)zd0 * (size_t)src_secsize + (size_t)yd0 * (size_t)nx1 + (size_t)xd0;
1820 float *dst = get_data() + (size_t)z0 * (size_t)dst_secsize + (size_t)y0 * (size_t)nx + (size_t)x0;
1821
1822 size_t src_gap = src_secsize - (y1-y0) * nx1;
1823 size_t dst_gap = dst_secsize - (y1-y0) * nx;
1824
1825 for (int i = z0; i < z1; i++) {
1826 for (int j = y0; j < y1; j++) {
1827 EMUtil::em_memcpy(dst, src, clipped_row_size);
1828 src += nx1;
1829 dst += nx;
1830 }
1831 src += src_gap;
1832 dst += dst_gap;
1833 }
1834
1835#ifdef EMAN2_USING_CUDA
1836 if(block->cudarwdata){
1837 copy_to_cuda(); // copy back to device as padding is faster on the host
1838 }
1839#endif
1840
1841 update();
1842 EXITFUNC;
1843}
1844
1845
1846void EMData::insert_scaled_sum(EMData *block, const FloatPoint &center,
1847 float scale, float mult_factor)
1848{
1849 ENTERFUNC;
1850 float * data = get_data();
1851 if (get_ndim()==3) {
1852 // Start by determining the region to operate on
1853 int xs=(int)floor(block->get_xsize()*scale/2.0);
1854 int ys=(int)floor(block->get_ysize()*scale/2.0);
1855 int zs=(int)floor(block->get_zsize()*scale/2.0);
1856 int x0=(int)center[0]-xs;
1857 int x1=(int)center[0]+xs;
1858 int y0=(int)center[1]-ys;
1859 int y1=(int)center[1]+ys;
1860 int z0=(int)center[2]-zs;
1861 int z1=(int)center[2]+zs;
1862
1863 if (x1<0||y1<0||z1<0||x0>get_xsize()||y0>get_ysize()||z0>get_zsize()) return; // object is completely outside the target volume
1864
1865 // make sure we stay inside the volume
1866 if (x0<0) x0=0;
1867 if (y0<0) y0=0;
1868 if (z0<0) z0=0;
1869 if (x1>=get_xsize()) x1=get_xsize()-1;
1870 if (y1>=get_ysize()) y1=get_ysize()-1;
1871 if (z1>=get_zsize()) z1=get_zsize()-1;
1872
1873 float bx=block->get_xsize()/2.0f;
1874 float by=block->get_ysize()/2.0f;
1875 float bz=block->get_zsize()/2.0f;
1876
1877 size_t idx;
1878 for (int x=x0; x<=x1; x++) {
1879 for (int y=y0; y<=y1; y++) {
1880 for (int z=z0; z<=z1; z++) {
1881 idx = x + y * nx + (size_t)z * nx * ny;
1882 if (scale==1) // skip interpolation so it is much faster
1883 data[idx] +=
1884 mult_factor*block->sget_value_at((x-center[0])+bx,(y-center[1])+by,(z-center[2])+bz);
1885 else
1886 data[idx] +=
1887 mult_factor*block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by,(z-center[2])/scale+bz);
1888
1889 }
1890 }
1891 }
1892 update();
1893 }
1894 else if (get_ndim()==2) {
1895 // Start by determining the region to operate on
1896 int xs=(int)floor(block->get_xsize()*scale/2.0);
1897 int ys=(int)floor(block->get_ysize()*scale/2.0);
1898 int x0=(int)center[0]-xs;
1899 int x1=(int)center[0]+xs;
1900 int y0=(int)center[1]-ys;
1901 int y1=(int)center[1]+ys;
1902
1903 if (x1<0||y1<0||x0>get_xsize()||y0>get_ysize()) return; // object is completely outside the target volume
1904
1905 // make sure we stay inside the volume
1906 if (x0<0) x0=0;
1907 if (y0<0) y0=0;
1908 if (x1>=get_xsize()) x1=get_xsize()-1;
1909 if (y1>=get_ysize()) y1=get_ysize()-1;
1910
1911 float bx=block->get_xsize()/2.0f;
1912 float by=block->get_ysize()/2.0f;
1913
1914 for (int x=x0; x<=x1; x++) {
1915 for (int y=y0; y<=y1; y++) {
1916 data[x + y * nx] +=
1917 mult_factor*block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by);
1918 }
1919 }
1920 update();
1921 }
1922 else {
1923 LOGERR("insert_scaled_sum supports only 2D and 3D data");
1924 throw ImageDimensionException("2D/3D only");
1925 }
1926
1927 EXITFUNC;
1928}
1929// else if ( m == 0 )
1930// {
1931// if ( n_f == -ny/2 )
1932// {
1933// t2++;
1934// // continue;
1935// for (int y = 0; y < return_slice->get_ysize(); ++y) {
1936// for (int x = 0; x < return_slice->get_xsize(); ++x) {
1937// double cur_val = return_slice->get_value_at(x,y);
1938// return_slice->set_value_at(x,y,cur_val+dat[idx]*std::pow(-1.0f,y));
1939// }
1940// }
1941// if (phase > 0.01 ) cout << "foo 2 " << phase << " " << amp << " " << dat[idx] << endl;
1942// }
1943// else
1944// {
1945// if ( n_f < 1 ) continue;
1946// t3++;
1947// for (int y = 0; y < return_slice->get_ysize(); ++y) {
1948// for (int x = 0; x < return_slice->get_xsize(); ++x) {
1949// double cur_val = return_slice->get_value_at(x,y);
1950// return_slice->set_value_at(x,y,cur_val+2*amp*cos(ndash*y+phase));
1951// }
1952// }
1953// }
1954// }
1955// else if ( n_f == -ny/2 )
1956// {
1957// if ( m == ((nx-2)/2) )
1958// {
1959// t4++;
1960// for (int y = 0; y < return_slice->get_ysize(); ++y) {
1961// for (int x = 0; x < return_slice->get_xsize(); ++x) {
1962// double cur_val = return_slice->get_value_at(x,y);
1963// return_slice->set_value_at(x,y,cur_val+dat[idx]*std::pow(-1.0f,x+y));
1964// }
1965// }
1966// if (phase > 0.01 ) cout << "foo 4 " << phase << " " << amp << " " << dat[idx] << endl;
1967// }
1968// else
1969// {
1970// t5++;
1971// for (int y = 0; y < return_slice->get_ysize(); ++y) {
1972// for (int x = 0; x < return_slice->get_xsize(); ++x) {
1973// double cur_val = return_slice->get_value_at(x,y);
1974// return_slice->set_value_at(x,y,cur_val+2*amp*cos(mdash*x+phase));
1975// }
1976// }
1977// }
1978// }
1979// else if ( n_f == 0 )
1980// {
1981// if ( m == ((nx-2)/2) )
1982// {
1983// t6++;
1984// for (int y = 0; y < return_slice->get_ysize(); ++y) {
1985// for (int x = 0; x < return_slice->get_xsize(); ++x) {
1986// double cur_val = return_slice->get_value_at(x,y);
1987// return_slice->set_value_at(x,y,cur_val+dat[idx]*std::pow(-1.0f,x));
1988// }
1989// }
1990// if (phase > 0.01 ) cout << "foo 3 " << phase << " " << amp << " " << dat[idx] << endl;
1991// }
1992// else
1993// {
1994// t7++;
1995// for (int y = 0; y < return_slice->get_ysize(); ++y) {
1996// for (int x = 0; x < return_slice->get_xsize(); ++x) {
1997// double cur_val = return_slice->get_value_at(x,y);
1998// return_slice->set_value_at(x,y,cur_val+2*amp*cos(mdash*x+M_PI*y + phase));
1999// }
2000// }
2001// }
2002// }
2003// else if ( m == ((nx-2)/2) )
2004// {
2005// if ( n_f < 1 ) continue;
2006// t8++;
2007// for (int y = 0; y < return_slice->get_ysize(); ++y) {
2008// for (int x = 0; x < return_slice->get_xsize(); ++x) {
2009// double cur_val = return_slice->get_value_at(x,y);
2010// return_slice->set_value_at(x,y,cur_val+2*amp*cos(ndash*y+M_PI*x+phase));
2011// }
2012// }
2013// }
2014// }
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::EMData::scale
void scale(float scale_factor)
scale the image by a factor.
Definition: emdata.cpp:886
EMAN::EMData::EMData
EMData()
For all image I/O.
Definition: emdata.cpp:70
EMAN::EMData::rdata
float * rdata
image real data
Definition: emdata.h:835
EMAN::EMData::attr_dict
Dict attr_dict
to store all image header info
Definition: emdata.h:833
EMAN::EMData::nx
int nx
image size
Definition: emdata.h:848
EMAN::EMData::flags
int flags
CTF data All CTF data become attribute ctf(vector<float>) in attr_dict –Grant Tang.
Definition: emdata.h:844
EMAN::EMData::nxyz
size_t nxyz
Definition: emdata.h:849
EMAN::EMUtil::em_memcpy
static void em_memcpy(void *dst, const void *const src, const size_t size)
Definition: emutil.h:384
EMAN::FloatPoint
FloatPoint defines a float-coordinate point in a 1D/2D/3D space.
Definition: geometry.h:278
EMAN::IntPoint
IntPoint defines an integer-coordinate point in a 1D/2D/3D space.
Definition: geometry.h:192
EMAN::Region
Region defines a 2D or 3D rectangular region specified by its origin coordinates and all edges' sizes...
Definition: geometry.h:497
EMAN::Util::ap2ri
static void ap2ri(float *data, size_t n)
convert complex data array from Amplitude/Phase format into Real/Imaginary format.
Definition: util.cpp:97
cuda_emfft.h
cuda_dd_fft_real_to_complex_nd
int cuda_dd_fft_real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz, int batch)
cuda_dd_fft_complex_to_real_nd
int cuda_dd_fft_complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz, int batch)
copy
EMData * copy() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
mult
void mult(int n)
multiply an integer number to each pixel value of the image.
Definition: emdata_core.h:97
set_value_at
void set_value_at(Vec3i loc, float val)
set_value_at with Vec3i
Definition: emdata_core.h:552
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
get_value_at
float get_value_at(int x, int y, int z) const
Get the pixel density value at coordinates (x,y,z).
Definition: emdata_core.h:221
copy_head
EMData * copy_head() const
Make an image with a copy of the current image's header.
Definition: emdata_core.cpp:100
write_image
void write_image(const string &filename, int img_index=0, EMUtil::ImageType imgtype=EMUtil::IMAGE_UNKNOWN, bool header_only=false, const Region *region=0, EMUtil::EMDataType filestoragetype=EMUtil::EM_FLOAT, bool use_host_endian=true)
write the header and data out to an image.
is_fftodd
bool is_fftodd() const
Does this image correspond to a (real-space) odd nx?
Definition: emdata_metadata.h:864
is_fftpadded
bool is_fftpadded() const
Is this image already extended along x for ffts?
Definition: emdata_metadata.h:831
set_fftodd
void set_fftodd(bool is_fftodd)
Mark this image as having (real-space) odd nx.
Definition: emdata_metadata.h:882
get_ysize
int get_ysize() const
Get the image y-dimensional size.
Definition: emdata_metadata.h:553
set_fftpad
void set_fftpad(bool is_fftpadded)
Mark this image as already extended along x for ffts.
Definition: emdata_metadata.h:850
set_ri
void set_ri(bool is_ri)
Mark this image as a real/imaginary format complex image.
Definition: emdata_metadata.h:817
is_ri
bool is_ri() const
Is this image a real/imaginary format complex image?
Definition: emdata_metadata.h:797
is_complex
bool is_complex() const
Is this a complex image?
Definition: emdata_metadata.h:648
get_zsize
int get_zsize() const
Get the image z-dimensional size.
Definition: emdata_metadata.h:562
get_xsize
int get_xsize() const
Get the image x-dimensional size.
Definition: emdata_metadata.h:544
get_ndim
int get_ndim() const
Get image dimension.
Definition: emdata_metadata.h:596
set_size
void set_size(int nx, int ny=1, int nz=1, bool noalloc=false)
Resize this EMData's main board memory pointer.
update
void update()
Mark EMData as changed, statistics, etc will be updated at need.
Definition: emdata_metadata.h:138
get_data
float * get_data() const
Get the image pixel density data in a 1D float array.
Definition: emdata_metadata.h:78
set_attr
void set_attr(const string &key, EMObject val)
Set a header attribute's value.
set_complex_x
void set_complex_x(bool is_complex_x)
Marks this image a 1D FFT image in X direction.
Definition: emdata_metadata.h:747
set_complex
void set_complex(bool is_complex)
Mark this image as a complex image.
Definition: emdata_metadata.h:712
vp
vector< pair< float, int > > vp
Definition: emdata_transform.cpp:53
calc_bessel
float calc_bessel(const int n, const float &x)
Definition: emdata_transform.cpp:1027
ri2inten
void ri2inten()
convert the complex image from real/imaginary to Intensity/0.
Definition: emdata_transform.cpp:979
do_fft_inplace
void do_fft_inplace()
Do FFT inplace.
Definition: emdata_transform.cpp:118
bispecRotTransInvDirect
EMData * bispecRotTransInvDirect(int type=0)
This computes the rotational and translational bispectral invariants of an image.
Definition: emdata_transform.cpp:1332
do_ift_inplace
void do_ift_inplace()
Definition: emdata_transform.cpp:392
render_ap24
EMBytes render_ap24(int x, int y, int xsize, int ysize, int bpl, float scale, int min_gray, int max_gray, float min_render, float max_render, float gamma, int flags)
Render the image into an 8-bit image.
Definition: emdata_transform.cpp:429
render_amp24
void render_amp24(int x, int y, int xsize, int ysize, int bpl, float scale, int min_gray, int max_gray, float min_render, float max_render, void *ref, void cmap(void *, int coord, unsigned char *tri))
Render the image into a 24-bit image.
Definition: emdata_transform.cpp:705
do_ift
EMData * do_ift()
return the inverse fourier transform (IFT) image of the current image.
Definition: emdata_transform.cpp:330
ap2ri
void ap2ri()
convert the complex image from amplitude/phase to real/imaginary
bispecRotTransInvN
EMData * bispecRotTransInvN(int N, int NK)
This computes the rotational and translational bispectral invariants of an image.
Definition: emdata_transform.cpp:1034
insert_scaled_sum
void insert_scaled_sum(EMData *block, const FloatPoint &center, float scale=1.0, float mult_factor=1.0)
Add a scaled image into another image at a specified location.
do_fft
EMData * do_fft() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
Definition: emdata_transform.cpp:55
ri2ap
void ri2ap()
convert the complex image from real/imaginary to amplitude/phase
Definition: emdata_transform.cpp:999
ImageFormatException
#define ImageFormatException(desc)
Definition: exception.h:147
ImageDimensionException
#define ImageDimensionException(desc)
Definition: exception.h:166
UnexpectedBehaviorException
#define UnexpectedBehaviorException(desc)
Definition: exception.h:400
insert_clip
void insert_clip(const EMData *const block, const IntPoint &origin)
Insert a clip into this image.
LOGWARN
#define LOGWARN
Definition: log.h:53
LOGERR
#define LOGERR
Definition: log.h:51
ENTERFUNC
#define ENTERFUNC
Definition: log.h:48
EXITFUNC
#define EXITFUNC
Definition: log.h:49
EMAN
E2Exception class.
Definition: aligner.h:40
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517
Generated by doxygen 1.9.2