EMAN2
emdata.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 "io/all_imageio.h"
34#include "ctf.h"
35#include "processor.h"
36#include "aligner.h"
37#include "cmp.h"
38#include "emfft.h"
39#include "projector.h"
40#include "geometry.h"
41#include <math.h>
42
43#include <gsl/gsl_sf_bessel.h>
44#include <gsl/gsl_errno.h>
45
46#include <iomanip>
47#include <complex>
48
49#include <algorithm> // fill
50#include <cmath>
51
52#ifdef WIN32
53 #define M_PI 3.14159265358979323846f
54#endif //WIN32
55
56#define EMDATA_EMAN2_DEBUG 0
57
58#ifdef EMAN2_USING_CUDA
59//#include <driver_functions.h>
60#include "cuda/cuda_processor.h"
61#include "cuda/cuda_emfft.h"
62#endif // EMAN2_USING_CUDA
63
64using namespace EMAN;
65using namespace std;
66using namespace boost;
67
68int EMData::totalalloc=0; // mainly used for debugging/memory leak purposes
69
70EMData::EMData() :
71#ifdef EMAN2_USING_CUDA
72 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
73#endif //EMAN2_USING_CUDA
74#ifdef FFT_CACHING
75 fftcache(0),
76#endif //FFT_CACHING
77 attr_dict(), rdata(0), supp(0), flags(0), changecount(0), nx(0), ny(0), nz(0), nxy(0), nxyz(0), xoff(0), yoff(0),
78 zoff(0), all_translation(), path(""), pathnum(0), rot_fp(0)
79
80{
81 ENTERFUNC;
82
83 attr_dict["apix_x"] = 1.0f;
84 attr_dict["apix_y"] = 1.0f;
85 attr_dict["apix_z"] = 1.0f;
86
87 attr_dict["is_complex"] = int(0);
88 attr_dict["is_complex_x"] = int(0);
89 attr_dict["is_complex_ri"] = int(1);
90
91 attr_dict["datatype"] = (int)EMUtil::EM_FLOAT;
92
93 EMData::totalalloc++;
94#ifdef MEMDEBUG2
95 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
96#endif
97
98 EXITFUNC;
99}
100
101EMData::EMData(const string& filename, int image_index) :
102#ifdef EMAN2_USING_CUDA
103 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
104#endif //EMAN2_USING_CUDA
105#ifdef FFT_CACHING
106 fftcache(0),
107#endif //FFT_CACHING
108 attr_dict(), rdata(0), supp(0), flags(0), changecount(0), nx(0), ny(0), nz(0), nxy(0), nxyz(0), xoff(0), yoff(0), zoff(0),
109 all_translation(), path(filename), pathnum(image_index), rot_fp(0)
110{
111 ENTERFUNC;
112
113 attr_dict["apix_x"] = 1.0f;
114 attr_dict["apix_y"] = 1.0f;
115 attr_dict["apix_z"] = 1.0f;
116
117 attr_dict["is_complex"] = int(0);
118 attr_dict["is_complex_x"] = int(0);
119 attr_dict["is_complex_ri"] = int(1);
120
121 attr_dict["datatype"] = (int)EMUtil::EM_FLOAT;
122
123 this->read_image(filename, image_index);
124
125 update();
126 EMData::totalalloc++;
127#ifdef MEMDEBUG2
128 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
129#endif
130
131 EXITFUNC;
132}
133
134EMData::EMData(const EMData& that) :
135#ifdef EMAN2_USING_CUDA
136 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
137#endif //EMAN2_USING_CUDA
138#ifdef FFT_CACHING
139 fftcache(0),
140#endif //FFT_CACHING
141 attr_dict(that.attr_dict), rdata(0), supp(0), flags(that.flags), changecount(that.changecount), nx(that.nx), ny(that.ny), nz(that.nz),
142 nxy(that.nx*that.ny), nxyz((size_t)that.nx*that.ny*that.nz), xoff(that.xoff), yoff(that.yoff), zoff(that.zoff),all_translation(that.all_translation), path(that.path),
143 pathnum(that.pathnum), rot_fp(0)
144{
145 ENTERFUNC;
146
147 float* data = that.rdata;
148 size_t num_bytes = (size_t)nx*ny*nz*sizeof(float);
149 if (data && num_bytes != 0)
150 {
151 rdata = (float*)EMUtil::em_malloc(num_bytes);
152 EMUtil::em_memcpy(rdata, data, num_bytes);
153 }
154#ifdef EMAN2_USING_CUDA
155 if (EMData::usecuda == 1 && num_bytes != 0 && that.cudarwdata != 0) {
156 //cout << "That copy constructor" << endl;
157 if(!rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");
158 cudaError_t error = cudaMemcpy(cudarwdata,that.cudarwdata,num_bytes,cudaMemcpyDeviceToDevice);
159 if ( error != cudaSuccess ) throw UnexpectedBehaviorException("cudaMemcpy failed in EMData copy construction with error: " + string(cudaGetErrorString(error)));
160 }
161#endif //EMAN2_USING_CUDA
162
163 if (that.rot_fp != 0) rot_fp = new EMData(*(that.rot_fp));
164
165 EMData::totalalloc++;
166#ifdef MEMDEBUG2
167 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
168#endif
169
170 ENTERFUNC;
171}
172
173EMData& EMData::operator=(const EMData& that)
174{
175 ENTERFUNC;
176
177 if ( this != &that )
178 {
179 free_memory(); // Free memory sets nx,ny and nz to 0
180
181 // Only copy the rdata if it exists, we could be in a scenario where only the header has been read
182 float* data = that.rdata;
183 size_t num_bytes = that.nx*that.ny*that.nz*sizeof(float);
184 if (data && num_bytes != 0)
185 {
186 nx = 1; // This prevents a memset in set_size
187 set_size(that.nx,that.ny,that.nz);
188 EMUtil::em_memcpy(rdata, data, num_bytes);
189 }
190
191 flags = that.flags;
192
193 all_translation = that.all_translation;
194
195 path = that.path;
196 pathnum = that.pathnum;
197 attr_dict = that.attr_dict;
198
199 xoff = that.xoff;
200 yoff = that.yoff;
201 zoff = that.zoff;
202
203#ifdef EMAN2_USING_CUDA
204 if (EMData::usecuda == 1 && num_bytes != 0 && that.cudarwdata != 0) {
205 if(cudarwdata){rw_free();}
206 if(!rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");;
207 cudaError_t error = cudaMemcpy(cudarwdata,that.cudarwdata,num_bytes,cudaMemcpyDeviceToDevice);
208 if ( error != cudaSuccess ) throw UnexpectedBehaviorException("cudaMemcpy failed in EMData copy construction with error: " + string(cudaGetErrorString(error)));
209 }
210#endif //EMAN2_USING_CUDA
211
212 changecount = that.changecount;
213
214 if (that.rot_fp != 0) rot_fp = new EMData(*(that.rot_fp));
215 else rot_fp = 0;
216 }
217 EXITFUNC;
218 return *this;
219}
220
221EMData::EMData(int nx, int ny, int nz, bool is_real) :
222#ifdef EMAN2_USING_CUDA
223 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
224#endif //EMAN2_USING_CUDA
225#ifdef FFT_CACHING
226 fftcache(0),
227#endif //FFT_CACHING
228 attr_dict(), rdata(0), supp(0), flags(0), changecount(0), nx(0), ny(0), nz(0), nxy(0), nxyz(0), xoff(0), yoff(0), zoff(0),
229 all_translation(), path(""), pathnum(0), rot_fp(0)
230{
231 ENTERFUNC;
232
233 // used to replace cube 'pixel'
234 attr_dict["apix_x"] = 1.0f;
235 attr_dict["apix_y"] = 1.0f;
236 attr_dict["apix_z"] = 1.0f;
237
238 if(is_real) { // create a real image [nx, ny, nz]
239 attr_dict["is_complex"] = int(0);
240 attr_dict["is_complex_x"] = int(0);
241 attr_dict["is_complex_ri"] = int(1);
242 set_size(nx, ny, nz);
243 }
244 else { //create a complex image which real dimension is [ny, ny, nz]
245 int new_nx = nx + 2 - nx%2;
246 set_size(new_nx, ny, nz);
247
248 attr_dict["is_complex"] = int(1);
249
250 if(ny==1 && nz ==1) {
251 attr_dict["is_complex_x"] = int(1);
252 }
253 else {
254 attr_dict["is_complex_x"] = int(0);
255 }
256
257 attr_dict["is_complex_ri"] = int(1);
258 attr_dict["is_fftpad"] = int(1);
259
260 if(nx%2 == 1) {
261 attr_dict["is_fftodd"] = 1;
262 }
263 }
264
265 this->to_zero();
266 update();
267 EMData::totalalloc++;
268#ifdef MEMDEBUG2
269 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
270#endif
271
272 EXITFUNC;
273}
274
275
276EMData::EMData(float* data, const int x, const int y, const int z, const Dict& attr_dict) :
277#ifdef EMAN2_USING_CUDA
278 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
279#endif //EMAN2_USING_CUDA
280#ifdef FFT_CACHING
281 fftcache(0),
282#endif //FFT_CACHING
283 attr_dict(attr_dict), rdata(data), supp(0), flags(0), changecount(0), nx(x), ny(y), nz(z), nxy(x*y), nxyz((size_t)x*y*z), xoff(0),
284 yoff(0), zoff(0), all_translation(), path(""), pathnum(0), rot_fp(0)
285{
286 ENTERFUNC;
287 // used to replace cube 'pixel'
288 attr_dict["apix_x"] = 1.0f;
289 attr_dict["apix_y"] = 1.0f;
290 attr_dict["apix_z"] = 1.0f;
291
292 EMData::totalalloc++;
293#ifdef MEMDEBUG2
294 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
295#endif
296
297 update();
298 EXITFUNC;
299}
300
301#ifdef EMAN2_USING_CUDA
302
303EMData::EMData(float* data, float* cudadata, const int x, const int y, const int z, const Dict& attr_dict) :
304 cudarwdata(cudadata), cudarodata(0), num_bytes(x*y*z*sizeof(float)), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(1),
305#ifdef FFT_CACHING
306 fftcache(0),
307#endif //FFT_CACHING
308 attr_dict(attr_dict), rdata(data), supp(0), flags(0), changecount(0), nx(x), ny(y), nz(z), nxy(x*y), nxyz((size_t)x*y*z), xoff(0),
309 yoff(0), zoff(0), all_translation(), path(""), pathnum(0), rot_fp(0)
310{
311 ENTERFUNC;
312
313 // used to replace cube 'pixel'
314 attr_dict["apix_x"] = 1.0f;
315 attr_dict["apix_y"] = 1.0f;
316 attr_dict["apix_z"] = 1.0f;
317
318 EMData::totalalloc++;
319#ifdef MEMDEBUG2
320 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
321#endif
322 update();
323 EXITFUNC;
324}
325
326#endif //EMAN2_USING_CUDA
327
328//debug
329using std::cout;
330using std::endl;
331EMData::~EMData()
332{
333 ENTERFUNC;
334#ifdef FFT_CACHING
335 if (fftcache!=0) { delete fftcache; fftcache=0;}
336#endif //FFT_CACHING
337 free_memory();
338
339#ifdef EMAN2_USING_CUDA
340 if(cudarwdata){rw_free();}
341 if(cudarodata){ro_free();}
342#endif // EMAN2_USING_CUDA
343 EMData::totalalloc--;
344#ifdef MEMDEBUG2
345 printf("EMDATA- %4d %p\n",EMData::totalalloc,this);
346#endif
347 EXITFUNC;
348}
349
350void EMData::clip_inplace(const Region & area,const float& fill_value)
351{
352 // Added by d.woolford
353 ENTERFUNC;
354
355// printf("cip %d %d %d %d %d %d %f %d %d %d\n",area.origin[0],area.origin[1],area.origin[2],area.size[0],area.size[1],area.size[2],fill_value,nx,ny,nz);
356 // Store the current dimension values
357 int prev_nx = nx, prev_ny = ny, prev_nz = nz;
358 size_t prev_size = (size_t)nx*ny*nz;
359
360 // Get the zsize, ysize and xsize of the final area, these are the new dimension sizes of the pixel data
361 int new_nz = ( area.size[2]==0 ? 1 : (int)area.size[2]);
362 int new_ny = ( area.size[1]==0 ? 1 : (int)area.size[1]);
363 int new_nx = (int)area.size[0];
364
365 if ( new_nz < 0 || new_ny < 0 || new_nx < 0 )
366 {
367 // Negative image dimensions were never tested nor considered when creating this implementation
368 throw ImageDimensionException("New image dimensions are negative - this is not supported in the clip_inplace operation");
369 }
370
371 size_t new_size = (size_t)new_nz*new_ny*new_nx;
372
373 // Get the translation values, they are used to construct the ClipInplaceVariables object
374 int x0 = (int) area.origin[0];
375 int y0 = (int) area.origin[1];
376 int z0 = (int) area.origin[2];
377
378 // Get a object that calculates all the interesting variables associated with the clip inplace operation
379 ClipInplaceVariables civ(prev_nx, prev_ny, prev_nz, new_nx, new_ny, new_nz, x0, y0, z0);
380
381 get_data(); // Do this here to make sure rdata is up to date, applicable if GPU stuff is occurring
382 // Double check to see if any memory shifting even has to occur
383 if ( x0 > prev_nx || y0 > prev_ny || z0 > prev_nz || civ.x_iter == 0 || civ.y_iter == 0 || civ.z_iter == 0)
384 {
385 // In this case the volume has been shifted beyond the location of the pixel rdata and
386 // the client should expect to see a volume with nothing in it.
387
388 // Set size calls realloc,
389 set_size(new_nx, new_ny, new_nz);
390
391 // Set pixel memory to zero - the client should expect to see nothing
392 EMUtil::em_memset(rdata, 0, (size_t)new_nx*new_ny*new_nz);
393
394 return;
395 }
396
397 // Resize the volume before memory shifting occurs if the new volume is larger than the previous volume
398 // All of the pixel rdata is guaranteed to be at the start of the new volume because realloc (called in set size)
399 // guarantees this.
400 if ( new_size > prev_size )
401 set_size(new_nx, new_ny, new_nz);
402
403 // Store the clipped row size.
404 size_t clipped_row_size = (civ.x_iter) * sizeof(float);
405
406 // Get the new sector sizes to save multiplication later.
407 size_t new_sec_size = new_nx * new_ny;
408 size_t prev_sec_size = prev_nx * prev_ny;
409
410 // Determine the memory locations of the source and destination pixels - at the point nearest
411 // to the beginning of the volume (rdata)
412 size_t src_it_begin = civ.prv_z_bottom*prev_sec_size + civ.prv_y_front*prev_nx + civ.prv_x_left;
413 size_t dst_it_begin = civ.new_z_bottom*new_sec_size + civ.new_y_front*new_nx + civ.new_x_left;
414
415 // This loop is in the forward direction (starting at points nearest to the beginning of the volume)
416 // it copies memory only when the destination pointer is less the source pointer - therefore
417 // ensuring that no memory "copied to" is yet to be "copied from"
418 for (int i = 0; i < civ.z_iter; ++i) {
419 for (int j = 0; j < civ.y_iter; ++j) {
420
421 // Determine the memory increments as dependent on i and j
422 // This could be optimized so that not so many multiplications are occurring...
423 size_t dst_inc = dst_it_begin + j*new_nx + i*new_sec_size;
424 size_t src_inc = src_it_begin + j*prev_nx + i*prev_sec_size;
425 float* local_dst = rdata + dst_inc;
426 float* local_src = rdata + src_inc;
427
428 if ( dst_inc >= src_inc )
429 {
430 // this is fine, it will happen now and then and it will be necessary to continue.
431 // the tempatation is to break, but you can't do that (because the point where memory intersects
432 // could be in this slice - and yes, this aspect could be optimized).
433 continue;
434 }
435
436 // Asserts are compiled only in debug mode
437 // This situation not encountered in testing thus far
438 Assert( dst_inc < new_size && src_inc < prev_size && dst_inc >= 0 && src_inc >= 0 );
439
440 // Finally copy the memory
441 memmove(local_dst, local_src, clipped_row_size);
442 }
443 }
444
445 // Determine the memory locations of the source and destination pixels - at the point nearest
446 // to the end of the volume (rdata+new_size)
447 size_t src_it_end = prev_size - civ.prv_z_top*prev_sec_size - civ.prv_y_back*prev_nx - prev_nx + civ.prv_x_left;
448 size_t dst_it_end = new_size - civ.new_z_top*new_sec_size - civ.new_y_back*new_nx - new_nx + civ.new_x_left;
449
450 // This loop is in the reverse direction (starting at points nearest to the end of the volume).
451 // It copies memory only when the destination pointer is greater than the source pointer therefore
452 // ensuring that no memory "copied to" is yet to be "copied from"
453 for (int i = 0; i < civ.z_iter; ++i) {
454 for (int j = 0; j < civ.y_iter; ++j) {
455
456 // Determine the memory increments as dependent on i and j
457 size_t dst_inc = dst_it_end - j*new_nx - i*new_sec_size;
458 size_t src_inc = src_it_end - j*prev_nx - i*prev_sec_size;
459 float* local_dst = rdata + dst_inc;
460 float* local_src = rdata + src_inc;
461
462 if (dst_inc <= (src_inc + civ.x_iter ))
463 {
464 // Overlap
465 if ( dst_inc > src_inc )
466 {
467 // Because the memcpy operation is the forward direction, and this "reverse
468 // direction" loop is proceeding in a backwards direction, it is possible
469 // that memory copied to is yet to be copied from (because memcpy goes forward).
470 // In this scenario pixel memory "copied to" is yet to be "copied from"
471 // i.e. there is overlap
472
473 // memmove handles overlapping cases.
474 // memmove could use a temporary buffer, or could go just go backwards
475 // the specification doesn't say how the function behaves...
476 // If memmove creates a temporary buffer is clip_inplace no longer inplace?
477 memmove(local_dst, local_src, clipped_row_size);
478 }
479 continue;
480 }
481
482 // This situation not encountered in testing thus far
483 Assert( dst_inc < new_size && src_inc < prev_size && dst_inc >= 0 && src_inc >= 0 );
484
485 // Perform the memory copy
486 EMUtil::em_memcpy(local_dst, local_src, clipped_row_size);
487 }
488 }
489
490 // Resize the volume after memory shifting occurs if the new volume is smaller than the previous volume
491 // set_size calls realloc, guaranteeing that the pixel rdata is in the right location.
492 if ( new_size < prev_size )
493 set_size(new_nx, new_ny, new_nz);
494
495 // Now set all the edges to zero
496
497 // Set the extra bottom z slices to the fill_value
498 if ( z0 < 0 )
499 {
500 //EMUtil::em_memset(rdata, 0, (-z0)*new_sec_size*sizeof(float));
501 size_t inc = (-z0)*new_sec_size;
502 std::fill(rdata,rdata+inc,fill_value);
503 }
504
505 // Set the extra top z slices to the fill_value
506 if ( civ.new_z_top > 0 )
507 {
508 float* begin_pointer = rdata + (new_nz-civ.new_z_top)*new_sec_size;
509 //EMUtil::em_memset(begin_pointer, 0, (civ.new_z_top)*new_sec_size*sizeof(float));
510 float* end_pointer = begin_pointer+(civ.new_z_top)*new_sec_size;
511 std::fill(begin_pointer,end_pointer,fill_value);
512 }
513
514 // Next deal with x and y edges by iterating through each slice
515 for ( int i = civ.new_z_bottom; i < civ.new_z_bottom + civ.z_iter; ++i )
516 {
517 // Set the extra front y components to the fill_value
518 if ( y0 < 0 )
519 {
520 float* begin_pointer = rdata + i*new_sec_size;
521 //EMUtil::em_memset(begin_pointer, 0, (-y0)*new_nx*sizeof(float));
522 float* end_pointer = begin_pointer+(-y0)*new_nx;
523 std::fill(begin_pointer,end_pointer,fill_value);
524 }
525
526 // Set the extra back y components to the fill_value
527 if ( civ.new_y_back > 0 )
528 {
529 float* begin_pointer = rdata + i*new_sec_size + (new_ny-civ.new_y_back)*new_nx;
530 //EMUtil::em_memset(begin_pointer, 0, (civ.new_y_back)*new_nx*sizeof(float));
531 float* end_pointer = begin_pointer+(civ.new_y_back)*new_nx;
532 std::fill(begin_pointer,end_pointer,fill_value);
533 }
534
535 // Iterate through the y to set each correct x component to the fill_value
536 for (int j = civ.new_y_front; j <civ.new_y_front + civ.y_iter; ++j)
537 {
538 // Set the extra left x components to the fill_value
539 if ( x0 < 0 )
540 {
541 float* begin_pointer = rdata + i*new_sec_size + j*new_nx;
542 //EMUtil::em_memset(begin_pointer, 0, (-x0)*sizeof(float));
543 float* end_pointer = begin_pointer+(-x0);
544 std::fill(begin_pointer,end_pointer,fill_value);
545 }
546
547 // Set the extra right x components to the fill_value
548 if ( civ.new_x_right > 0 )
549 {
550 float* begin_pointer = rdata + i*new_sec_size + j*new_nx + (new_nx - civ.new_x_right);
551 //EMUtil::em_memset(begin_pointer, 0, (civ.new_x_right)*sizeof(float));
552 float* end_pointer = begin_pointer+(civ.new_x_right);
553 std::fill(begin_pointer,end_pointer,fill_value);
554 }
555
556 }
557 }
558
559// These couts may be useful
560// cout << "start starts " << civ.prv_x_left << " " << civ.prv_y_front << " " << civ.prv_z_bottom << endl;
561// cout << "start ends " << civ.prv_x_right << " " << civ.prv_y_back << " " << civ.prv_z_top << endl;
562// cout << "dst starts " << civ.new_x_left << " " << civ.new_y_front << " " << civ.new_z_bottom << endl;
563// cout << "dst ends " << civ.new_x_right << " " << civ.new_y_back << " " << civ.new_z_top << endl;
564// cout << "total iter z - " << civ.z_iter << " y - " << civ.y_iter << " x - " << civ.x_iter << endl;
565// cout << "=====" << endl;
566// cout << "dst_end is " << dst_it_end << " src end is " << src_it_end << endl;
567// cout << "dst_begin is " << dst_it_begin << " src begin is " << src_it_begin << endl;
568
569 // Update appropriate attributes (Copied and pasted from get_clip)
570 if( attr_dict.has_key("origin_x") && attr_dict.has_key("origin_y") &&
571 attr_dict.has_key("origin_z") )
572 {
573 float xorigin = attr_dict["origin_x"];
574 float yorigin = attr_dict["origin_y"];
575 float zorigin = attr_dict["origin_z"];
576
577 float apix_x = attr_dict["apix_x"];
578 float apix_y = attr_dict["apix_y"];
579 float apix_z = attr_dict["apix_z"];
580
581 set_xyz_origin(xorigin + apix_x * area.origin[0],
582 yorigin + apix_y * area.origin[1],
583 zorigin + apix_z * area.origin[2]);
584 }
585
586 // Set the update flag because the size of the image has changed and stats should probably be recalculated if requested.
587 update();
588
589 EXITFUNC;
590}
591
592EMData *EMData::get_clip(const Region & area, const float fill) const
593{
594 ENTERFUNC;
595 // Steve: removed check 10/1/17. Prevented some sane operations, like a 3d clip which was actually 2d from a 2d volume
596// if (get_ndim() != area.get_ndim()) {
597// LOGERR("cannot get %dD clip out of %dD image", area.get_ndim(),get_ndim());
598// return 0;
599// }
600
601 EMData *result = new EMData();
602
603 // Ensure that all of the metadata of this is stored in the new object
604 // Originally added to ensure that euler angles were retained when preprocessing (zero padding) images
605 // prior to insertion into the 3D for volume in the reconstruction phase (see reconstructor.cpp/h).
606 result->attr_dict = this->attr_dict;
607 int zsize = (int)area.size[2];
608 if (zsize == 0 && nz <= 1) {
609 zsize = 1;
610 }
611 int ysize = (ny<=1 && (int)area.size[1]==0 ? 1 : (int)area.size[1]);
612
613 if ( (int)area.size[0] < 0 || ysize < 0 || zsize < 0 )
614 {
615 // Negative image dimensions not supported - added retrospectively by d.woolford (who didn't write get_clip but wrote clip_inplace)
616 throw ImageDimensionException("New image dimensions are negative - this is not supported in the the get_clip operation");
617 }
618
619//#ifdef EMAN2_USING_CUDA
620 // Strategy is always to prefer using the GPU if possible
621// bool use_gpu = false;
622// if ( gpu_operation_preferred() ) {
623// result->set_size_cuda((int)area.size[0], ysize, zsize);
624 //CudaDataLock lock(this); // Just so we never have to recopy this data to and from the GPU
625// result->get_cuda_data(); // Force the allocation - set_size_cuda is lazy
626 // Setting the value is necessary seeing as cuda data is not automatically zeroed
627// result->to_value(fill); // This will automatically use the GPU.
628// use_gpu = true;
629// } else { // cpu == True
630// result->set_size((int)area.size[0], ysize, zsize);
631// if (fill != 0.0) { result->to_value(fill); };
632// }
633//#else
634 result->set_size((int)area.size[0], ysize, zsize);
635 if (fill != 0.0) { result->to_value(fill); };
636//#endif //EMAN2_USING_CUDA
637
638 int x0 = (int) area.origin[0];
639 x0 = x0 < 0 ? 0 : x0;
640
641 int y0 = (int) area.origin[1];
642 y0 = y0 < 0 ? 0 : y0;
643
644 int z0 = (int) area.origin[2];
645 z0 = z0 < 0 ? 0 : z0;
646
647 int x1 = (int) (area.origin[0] + area.size[0]);
648 x1 = x1 > nx ? nx : x1;
649
650 int y1 = (int) (area.origin[1] + area.size[1]);
651 y1 = y1 > ny ? ny : y1;
652
653 int z1 = (int) (area.origin[2] + area.size[2]);
654 z1 = z1 > nz ? nz : z1;
655 if (z1 <= 0) {
656 z1 = 1;
657 }
658
659 result->insert_clip(this,-((IntPoint)area.origin));
660
661 if( attr_dict.has_key("apix_x") && attr_dict.has_key("apix_y") &&
662 attr_dict.has_key("apix_z") )
663 {
664 float apix_x = attr_dict["apix_x"];
665 float apix_y = attr_dict["apix_y"];
666 float apix_z = attr_dict["apix_z"];
667 result->set_attr("apix_x",apix_x);
668 result->set_attr("apix_y",apix_y);
669 result->set_attr("apix_z",apix_z);
670
671 if( attr_dict.has_key("origin_x") && attr_dict.has_key("origin_y") &&
672 attr_dict.has_key("origin_z") )
673 {
674 float xorigin = attr_dict["origin_x"];
675 float yorigin = attr_dict["origin_y"];
676 float zorigin = attr_dict["origin_z"];
677
678
679 result->set_xyz_origin(xorigin + apix_x * area.origin[0],
680 yorigin + apix_y * area.origin[1],
681 zorigin + apix_z * area.origin[2]);
682 }
683 else {
684 result->set_xyz_origin(apix_x * area.origin[0],
685 apix_y * area.origin[1],
686 apix_z * area.origin[2]);
687 }
688
689 }
690
691//#ifdef EMAN2_USING_CUDA
692// if (use_gpu) result->gpu_update();
693// else result->update();
694//#else
695 result->update();
696//#endif // EMAN2_USING_CUDA
697
698
699 result->set_path(path);
700 result->set_pathnum(pathnum);
701
702 EXITFUNC;
703 return result;
704}
705
706
707EMData *EMData::get_top_half() const
708{
709 ENTERFUNC;
710
711 if (get_ndim() != 3) {
712 throw ImageDimensionException("3D only");
713 }
714
715 EMData *half = new EMData();
716 half->attr_dict = attr_dict;
717 half->set_size(nx, ny, nz / 2);
718
719 float *half_data = half->get_data();
720 EMUtil::em_memcpy(half_data, &(get_data()[(size_t)nz / 2 * (size_t)nx * (size_t)ny]), sizeof(float) * (size_t)nx * (size_t)ny * (size_t)nz / 2lu);
721
722 float apix_z = attr_dict["apix_z"];
723 float origin_z = attr_dict["origin_z"];
724 origin_z += apix_z * nz / 2;
725 half->attr_dict["origin_z"] = origin_z;
726 half->update();
727
728 EXITFUNC;
729 return half;
730}
731
732
733EMData *EMData::get_rotated_clip(const Transform &xform,
734 const IntSize &size, float)
735{
736 EMData *result = new EMData();
737 result->set_size(size[0],size[1],size[2]);
738
739 if (nz==1) {
740 for (int y=-size[1]/2; y<(size[1]+1)/2; y++) {
741 for (int x=-size[0]/2; x<(size[0]+1)/2; x++) {
742 Vec3f xv=xform.transform(Vec3f((float)x,(float)y,0.0f));
743 float v = 0;
744
745 if (xv[0]<0||xv[1]<0||xv[0]>nx-2||xv[1]>ny-2) v=0.;
746 else v=sget_value_at_interp(xv[0],xv[1]);
747 result->set_value_at(x+size[0]/2,y+size[1]/2,v);
748 }
749 }
750 }
751 else {
752 for (int z=-size[2]/2; z<(size[2]+1)/2; z++) {
753 for (int y=-size[1]/2; y<(size[1]+1)/2; y++) {
754 for (int x=-size[0]/2; x<(size[0]+1)/2; x++) {
755 Vec3f xv=xform.transform(Vec3f((float)x,(float)y,(float)z));
756 float v = 0;
757
758 if (xv[0]<0||xv[1]<0||xv[2]<0||xv[0]>nx-2||xv[1]>ny-2||xv[2]>nz-2) v=0.;
759 else v=sget_value_at_interp(xv[0],xv[1],xv[2]);
760 result->set_value_at(x+size[0]/2,y+size[1]/2,z+size[2]/2,v);
761 }
762 }
763 }
764 }
765
766 result->attr_dict["apix_x"] = attr_dict["apix_x"];
767 result->attr_dict["apix_y"] = attr_dict["apix_y"];
768 result->attr_dict["apix_z"] = attr_dict["apix_z"];
769
770 result->update();
771
772 return result;
773}
774
775
776EMData* EMData::window_center(int l) {
777 ENTERFUNC;
778 // sanity checks
779 int n = nx;
780 if (is_complex()) {
781 LOGERR("Need real-space data for window_center()");
782 throw ImageFormatException(
783 "Complex input image; real-space expected.");
784 }
785 if (is_fftpadded()) {
786 // image has been fft-padded, compute the real-space size
787 n -= (2 - int(is_fftodd()));
788 }
789 int corner = n/2 - l/2;
790 int ndim = get_ndim();
791 EMData* ret;
792 switch (ndim) {
793 case 3:
794 if ((n != ny) || (n != nz)) {
795 LOGERR("Need the real-space image to be cubic.");
796 throw ImageFormatException(
797 "Need cubic real-space image.");
798 }
799 ret = get_clip(Region(corner, corner, corner, l, l, l));
800 break;
801 case 2:
802 if (n != ny) {
803 LOGERR("Need the real-space image to be square.");
804 throw ImageFormatException(
805 "Need square real-space image.");
806 }
807 //cout << "Using corner " << corner << endl;
808 ret = get_clip(Region(corner, corner, l, l));
809 break;
810 case 1:
811 ret = get_clip(Region(corner, l));
812 break;
813 default:
814 throw ImageDimensionException(
815 "window_center only supports 1-d, 2-d, and 3-d images");
816 }
817 return ret;
818 EXITFUNC;
819}
820
821
822float *EMData::setup4slice(bool redo)
823{
824 ENTERFUNC;
825
826 if (!is_complex()) {
827 throw ImageFormatException("complex image only");
828 }
829
830 if (get_ndim() != 3) {
831 throw ImageDimensionException("3D only");
832 }
833
834 if (supp) {
835 if (redo) {
836 EMUtil::em_free(supp);
837 supp = 0;
838 }
839 else {
840 EXITFUNC;
841 return supp;
842 }
843 }
844
845 const int SUPP_ROW_SIZE = 8;
846 const int SUPP_ROW_OFFSET = 4;
847 const int supp_size = SUPP_ROW_SIZE + SUPP_ROW_OFFSET;
848
849 supp = (float *) EMUtil::em_calloc(supp_size * ny * nz, sizeof(float));
850 int nxy = nx * ny;
851 int supp_xy = supp_size * ny;
852 float * data = get_data();
853
854 for (int z = 0; z < nz; z++) {
855 size_t cur_z1 = z * nxy;
856 size_t cur_z2 = z * supp_xy;
857
858 for (int y = 0; y < ny; y++) {
859 size_t cur_y1 = y * nx;
860 size_t cur_y2 = y * supp_size;
861
862 for (int x = 0; x < SUPP_ROW_SIZE; x++) {
863 size_t k = (x + SUPP_ROW_OFFSET) + cur_y2 + cur_z2;
864 supp[k] = data[x + cur_y1 + cur_z1];
865 }
866 }
867 }
868
869 for (int z = 1, zz = nz - 1; z < nz; z++, zz--) {
870 size_t cur_z1 = zz * nxy;
871 size_t cur_z2 = z * supp_xy;
872
873 for (int y = 1, yy = ny - 1; y < ny; y++, yy--) {
874 supp[y * 12 + cur_z2] = data[4 + yy * nx + cur_z1];
875 supp[1 + y * 12 + cur_z2] = -data[5 + yy * nx + cur_z1];
876 supp[2 + y * 12 + cur_z2] = data[2 + yy * nx + cur_z1];
877 supp[3 + y * 12 + cur_z2] = -data[3 + yy * nx + cur_z1];
878 }
879 }
880
881 EXITFUNC;
882 return supp;
883}
884
885
886void EMData::scale(float s)
887{
888 ENTERFUNC;
889 Transform t;
890 t.set_scale(s);
891 transform(t);
892 EXITFUNC;
893}
894
895
896void EMData::translate(int dx, int dy, int dz)
897{
898 ENTERFUNC;
899 translate(Vec3i(dx, dy, dz));
900 EXITFUNC;
901}
902
903
904void EMData::translate(float dx, float dy, float dz)
905{
906 ENTERFUNC;
907 int dx_ = Util::round(dx);
908 int dy_ = Util::round(dy);
909 int dz_ = Util::round(dz);
910 if( ( (dx-dx_) == 0 ) && ( (dy-dy_) == 0 ) && ( (dz-dz_) == 0 )) {
911 translate(dx_, dy_, dz_);
912 }
913 else {
914 translate(Vec3f(dx, dy, dz));
915 }
916 EXITFUNC;
917}
918
919
920void EMData::translate(const Vec3i &translation)
921{
922 ENTERFUNC;
923
924 //if traslation is 0, do nothing
925 if( translation[0] == 0 && translation[1] == 0 && translation[2] == 0) {
926 EXITFUNC;
927 return;
928 }
929
930 Dict params("trans",static_cast< vector<int> >(translation));
931 process_inplace("xform.translate.int",params);
932
933 // update() - clip_inplace does the update
934 all_translation += translation;
935
936 EXITFUNC;
937}
938
939
940void EMData::translate(const Vec3f &translation)
941{
942 ENTERFUNC;
943
944 if( translation[0] == 0.0f && translation[1] == 0.0f && translation[2] == 0.0f ) {
945 EXITFUNC;
946 return;
947 }
948
949 Transform* t = new Transform();
950 t->set_trans(translation);
951 process_inplace("xform",Dict("transform",t));
952 delete t;
953
954 all_translation += translation;
955 EXITFUNC;
956}
957
958
959void EMData::rotate(float az, float alt, float phi)
960{
961 Dict d("type","eman");
962 d["az"] = az;
963 d["alt"] = alt;
964 d["phi"] = phi;
965 Transform t(d);
966 transform(t);
967}
968
969
970
971void EMData::rotate(const Transform & t)
972{
973 cout << "Deprecation warning in EMData::rotate. Please consider using EMData::transform() instead " << endl;
974 transform(t);
975}
976
977float EMData::max_3D_pixel_error(const Transform &t1, const Transform & t2, float r) {
978
979 Transform t;
980 int r0 = (int)r;
981 float ddmax = 0.0f;
982
983 t = t2*t1.inverse();
984
985 for (int i=0; i<int(2*M_PI*r0+0.5); i++) {
986 float ang = (float)i/r;
987 Vec3f v = Vec3f(r0*cos(ang), r0*sin(ang), 0.0f);
988 Vec3f d = t*v-v;
989 ddmax = Util::get_max(ddmax, d[0]*d[0]+d[1]*d[1]+d[2]*d[2]);
990 }
991
992 return std::sqrt(ddmax);
993}
994
995void EMData::rotate_translate(float az, float alt, float phi, float dx, float dy, float dz)
996{
997 cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
998// Transform3D t( az, alt, phi,Vec3f(dx, dy, dz));
999 Transform t;
1000 t.set_rotation(Dict("type", "eman", "az", az, "alt", alt, "phi", phi));
1001 t.set_trans(dx, dy, dz);
1002 rotate_translate(t);
1003}
1004
1005
1006void EMData::rotate_translate(float az, float alt, float phi, float dx, float dy,
1007 float dz, float pdx, float pdy, float pdz)
1008{
1009 cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
1010// Transform3D t(Vec3f(dx, dy, dz), az, alt, phi, Vec3f(pdx,pdy,pdz));
1011// rotate_translate(t);
1012
1013 Transform t;
1014 t.set_pre_trans(Vec3f(dx, dy, dz));
1015 t.set_rotation(Dict("type", "eman", "az", az, "alt", alt, "phi", phi));
1016 t.set_trans(pdx, pdy, pdz);
1017 rotate_translate(t);
1018}
1019
1020//void EMData::rotate_translate(const Transform3D & RA)
1021//{
1022// cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
1023// ENTERFUNC;
1024//
1025//#if EMDATA_EMAN2_DEBUG
1026// std::cout << "start rotate_translate..." << std::endl;
1027//#endif
1028//
1029// float scale = RA.get_scale();
1030// Vec3f dcenter = RA.get_center();
1031// Vec3f translation = RA.get_posttrans();
1032// Dict rotation = RA.get_rotation(Transform3D::EMAN);
1034// Transform3D RAInv = RA.inverse(); // We're rotating the coordinate system, not the data
1036//#if EMDATA_EMAN2_DEBUG
1037// vector<string> keys = rotation.keys();
1038// vector<string>::const_iterator it;
1039// for(it=keys.begin(); it!=keys.end(); ++it) {
1041// std::cout << *it << " : " << (float)rotation.get(*it) << std::endl;
1042// }
1043//#endif
1044// float inv_scale = 1.0f;
1045//
1046// if (scale != 0) {
1047// inv_scale = 1.0f / scale;
1048// }
1049//
1050// float *src_data = 0;
1051// float *des_data = 0;
1052//
1053// src_data = get_data();
1054// des_data = (float *) EMUtil::em_malloc(nx * ny * nz * sizeof(float));
1055//
1056// if (nz == 1) {
1057// float x2c = nx / 2 - dcenter[0] + RAInv[0][3];
1058// float y2c = ny / 2 - dcenter[1] + RAInv[1][3];
1059// float y = -ny / 2 + dcenter[1]; // changed 0 to 1 in dcenter and below
1060// for (int j = 0; j < ny; j++, y += 1.0f) {
1061// float x = -nx / 2 + dcenter[0];
1062// for (int i = 0; i < nx; i++, x += 1.0f) {
1063// float x2 = RAInv[0][0]*x + RAInv[0][1]*y + x2c;
1064// float y2 = RAInv[1][0]*x + RAInv[1][1]*y + y2c;
1065//
1066// if (x2 < 0 || x2 >= nx || y2 < 0 || y2 >= ny ) {
1067// des_data[i + j * nx] = 0; // It may be tempting to set this value to the
1068// // mean but in fact this is not a good thing to do. Talk to S.Ludtke about it.
1069// }
1070// else {
1071// int ii = Util::fast_floor(x2);
1072// int jj = Util::fast_floor(y2);
1073// int k0 = ii + jj * nx;
1074// int k1 = k0 + 1;
1075// int k2 = k0 + nx;
1076// int k3 = k0 + nx + 1;
1077//
1078// if (ii == nx - 1) {
1079// k1--;
1080// k3--;
1081// }
1082// if (jj == ny - 1) {
1083// k2 -= nx;
1084// k3 -= nx;
1085// }
1086//
1087// float t = x2 - ii;
1088// float u = y2 - jj;
1089//
1090// des_data[i + j * nx] = Util::bilinear_interpolate(src_data[k0],src_data[k1], src_data[k2], src_data[k3],t,u); // This is essentially basic interpolation
1091// }
1092// }
1093// }
1094// }
1095// else {
1096//#if EMDATA_EMAN2_DEBUG
1097// std::cout << "This is the 3D case." << std::endl ;
1098//#endif
1099//
1100// Transform3D mx = RA;
1101// mx.set_scale(inv_scale);
1102//
1103//#if EMDATA_EMAN2_DEBUG
1106//#endif
1107//
1108// int nxy = nx * ny;
1109// int l = 0;
1110//
1111// float x2c = nx / 2 - dcenter[0] + RAInv[0][3];;
1112// float y2c = ny / 2 - dcenter[1] + RAInv[1][3];;
1113// float z2c = nz / 2 - dcenter[2] + RAInv[2][3];;
1114//
1115// float z = -nz / 2 + dcenter[2]; //
1116//
1117// size_t ii, k0, k1, k2, k3, k4, k5, k6, k7;
1118// for (int k = 0; k < nz; k++, z += 1.0f) {
1119// float y = -ny / 2 + dcenter[1]; //
1120// for (int j = 0; j < ny; j++, y += 1.0f) {
1121// float x = -nx / 2 + dcenter[0];
1122// for (int i = 0; i < nx; i++, l++ , x += 1.0f) {
1123// float x2 = RAInv[0][0] * x + RAInv[0][1] * y + RAInv[0][2] * z + x2c;
1124// float y2 = RAInv[1][0] * x + RAInv[1][1] * y + RAInv[1][2] * z + y2c;
1125// float z2 = RAInv[2][0] * x + RAInv[2][1] * y + RAInv[2][2] * z + z2c;
1126//
1127//
1128// if (x2 < 0 || y2 < 0 || z2 < 0 ||
1129// x2 >= nx || y2 >= ny || z2>= nz ) {
1130// des_data[l] = 0;
1131// }
1132// else {
1133// int ix = Util::fast_floor(x2);
1134// int iy = Util::fast_floor(y2);
1135// int iz = Util::fast_floor(z2);
1136// float tuvx = x2-ix;
1137// float tuvy = y2-iy;
1138// float tuvz = z2-iz;
1139// ii = ix + iy * nx + iz * nxy;
1140//
1141// k0 = ii;
1142// k1 = k0 + 1;
1143// k2 = k0 + nx;
1144// k3 = k0 + nx+1;
1145// k4 = k0 + nxy;
1146// k5 = k1 + nxy;
1147// k6 = k2 + nxy;
1148// k7 = k3 + nxy;
1149//
1150// if (ix == nx - 1) {
1151// k1--;
1152// k3--;
1153// k5--;
1154// k7--;
1155// }
1156// if (iy == ny - 1) {
1157// k2 -= nx;
1158// k3 -= nx;
1159// k6 -= nx;
1160// k7 -= nx;
1161// }
1162// if (iz == nz - 1) {
1163// k4 -= nxy;
1164// k5 -= nxy;
1165// k6 -= nxy;
1166// k7 -= nxy;
1167// }
1168//
1169// des_data[l] = Util::trilinear_interpolate(src_data[k0],
1170// src_data[k1],
1171// src_data[k2],
1172// src_data[k3],
1173// src_data[k4],
1174// src_data[k5],
1175// src_data[k6],
1176// src_data[k7],
1177// tuvx, tuvy, tuvz);
1178//#if EMDATA_EMAN2_DEBUG
1179// printf(" ix=%d \t iy=%d \t iz=%d \t value=%f \n", ix ,iy, iz, des_data[l] );
1180// std::cout << src_data[ii] << std::endl;
1181//#endif
1182// }
1183// }
1184// }
1185// }
1186// }
1187//
1188// if( rdata )
1189// {
1190// EMUtil::em_free(rdata);
1191// rdata = 0;
1192// }
1193// rdata = des_data;
1194//
1195// scale_pixel(inv_scale);
1196//
1197// attr_dict["origin_x"] = (float) attr_dict["origin_x"] * inv_scale;
1198// attr_dict["origin_y"] = (float) attr_dict["origin_y"] * inv_scale;
1199// attr_dict["origin_z"] = (float) attr_dict["origin_z"] * inv_scale;
1200//
1201// update();
1202// all_translation += translation;
1203// EXITFUNC;
1204//}
1205
1206
1207
1208
1209void EMData::rotate_x(int dx)
1210{
1211 ENTERFUNC;
1212
1213 if (get_ndim() > 2) {
1214 throw ImageDimensionException("no 3D image");
1215 }
1216
1217
1218 size_t row_size = nx * sizeof(float);
1219 float *tmp = (float*)EMUtil::em_malloc(row_size);
1220 float * data = get_data();
1221
1222 for (int y = 0; y < ny; y++) {
1223 int y_nx = y * nx;
1224 for (int x = 0; x < nx; x++) {
1225 tmp[x] = data[y_nx + (x + dx) % nx];
1226 }
1227 EMUtil::em_memcpy(&data[y_nx], tmp, row_size);
1228 }
1229
1230 update();
1231 if( tmp )
1232 {
1233 delete[]tmp;
1234 tmp = 0;
1235 }
1236 EXITFUNC;
1237}
1238
1239double EMData::dot_rotate_translate(EMData * with, float dx, float dy, float da, const bool mirror)
1240{
1241 ENTERFUNC;
1242
1243 if (!EMUtil::is_same_size(this, with)) {
1244 LOGERR("images not same size");
1245 throw ImageFormatException("images not same size");
1246 }
1247
1248 if (get_ndim() == 3) {
1249 LOGERR("1D/2D Images only");
1250 throw ImageDimensionException("1D/2D only");
1251 }
1252
1253 float *this_data = 0;
1254
1255 this_data = get_data();
1256
1257 float da_rad = da*(float)M_PI/180.0f;
1258
1259 float *with_data = with->get_data();
1260 float mx0 = cos(da_rad);
1261 float mx1 = sin(da_rad);
1262 float y = -ny / 2.0f;
1263 float my0 = mx0 * (-nx / 2.0f - 1.0f) + nx / 2.0f - dx;
1264 float my1 = -mx1 * (-nx / 2.0f - 1.0f) + ny / 2.0f - dy;
1265 double result = 0;
1266
1267 for (int j = 0; j < ny; j++) {
1268 float x2 = my0 + mx1 * y;
1269 float y2 = my1 + mx0 * y;
1270
1271 int ii = Util::fast_floor(x2);
1272 int jj = Util::fast_floor(y2);
1273 float t = x2 - ii;
1274 float u = y2 - jj;
1275
1276 for (int i = 0; i < nx; i++) {
1277 t += mx0;
1278 u -= mx1;
1279
1280 if (t >= 1.0f) {
1281 ii++;
1282 t -= 1.0f;
1283 }
1284
1285 if (u >= 1.0f) {
1286 jj++;
1287 u -= 1.0f;
1288 }
1289
1290 if (t < 0) {
1291 ii--;
1292 t += 1.0f;
1293 }
1294
1295 if (u < 0) {
1296 jj--;
1297 u += 1.0f;
1298 }
1299
1300 if (ii >= 0 && ii <= nx - 2 && jj >= 0 && jj <= ny - 2) {
1301 int k0 = ii + jj * nx;
1302 int k1 = k0 + 1;
1303 int k2 = k0 + nx + 1;
1304 int k3 = k0 + nx;
1305
1306 float tt = 1 - t;
1307 float uu = 1 - u;
1308 int idx = i + j * nx;
1309 if (mirror) idx = nx-1-i+j*nx; // mirroring of Transforms is always about the y axis
1310 result += (this_data[k0] * tt * uu + this_data[k1] * t * uu +
1311 this_data[k2] * t * u + this_data[k3] * tt * u) * with_data[idx];
1312 }
1313 }
1314 y += 1.0f;
1315 }
1316
1317 EXITFUNC;
1318 return result;
1319}
1320
1321
1322EMData *EMData::little_big_dot(EMData * with, bool do_sigma)
1323{
1324 ENTERFUNC;
1325
1326 if (get_ndim() > 2) {
1327 throw ImageDimensionException("1D/2D only");
1328 }
1329
1330 EMData *ret = copy_head();
1331 ret->set_size(nx,ny,nz);
1332 ret->to_zero();
1333
1334 int nx2 = with->get_xsize();
1335 int ny2 = with->get_ysize();
1336 float em = with->get_edge_mean();
1337
1338 float *data = get_data();
1339 float *with_data = with->get_data();
1340 float *ret_data = ret->get_data();
1341
1342 float sum2 = (Util::square((float)with->get_attr("sigma")) +
1343 Util::square((float)with->get_attr("mean")));
1344 if (do_sigma) {
1345 for (int j = ny2 / 2; j < ny - ny2 / 2; j++) {
1346 for (int i = nx2 / 2; i < nx - nx2 / 2; i++) {
1347 float sum = 0;
1348 float sum1 = 0;
1349 float summ = 0;
1350 int k = 0;
1351
1352 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
1353 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
1354 int l = ii + jj * nx;
1355 sum1 += Util::square(data[l]);
1356 summ += data[l];
1357 sum += data[l] * with_data[k];
1358 k++;
1359 }
1360 }
1361 float tmp_f1 = (sum1 / 2.0f - sum) / (nx2 * ny2);
1362 float tmp_f2 = Util::square((float)with->get_attr("mean") -
1363 summ / (nx2 * ny2));
1364 ret_data[i + j * nx] = sum2 + tmp_f1 - tmp_f2;
1365 }
1366 }
1367 }
1368 else {
1369 for (int j = ny2 / 2; j < ny - ny2 / 2; j++) {
1370 for (int i = nx2 / 2; i < nx - nx2 / 2; i++) {
1371 float eml = 0;
1372 float dot = 0;
1373 float dot2 = 0;
1374
1375 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
1376 eml += data[ii + (j - ny2 / 2) * nx] + data[ii + (j + ny2 / 2 - 1) * nx];
1377 }
1378
1379 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
1380 eml += data[i - nx2 / 2 + jj * nx] + data[i + nx2 / 2 - 1 + jj * nx];
1381 }
1382
1383 eml /= (nx2 + ny2) * 2.0f;
1384 int k = 0;
1385
1386 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
1387 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
1388 dot += (data[ii + jj * nx] - eml) * (with_data[k] - em);
1389 dot2 += Util::square(data[ii + jj * nx] - eml);
1390 k++;
1391 }
1392 }
1393
1394 dot2 = std::sqrt(dot2);
1395
1396 if (dot2 == 0) {
1397 ret_data[i + j * nx] = 0;
1398 }
1399 else {
1400 ret_data[i + j * nx] = dot / (nx2 * ny2 * dot2 * (float)with->get_attr("sigma"));
1401 }
1402 }
1403 }
1404 }
1405
1406 ret->update();
1407
1408 EXITFUNC;
1409 return ret;
1410}
1411
1412
1413EMData *EMData::do_radon()
1414{
1415 ENTERFUNC;
1416
1417 if (get_ndim() != 2) {
1418 throw ImageDimensionException("2D only");
1419 }
1420
1421 if (nx != ny) {
1422 throw ImageFormatException("square image only");
1423 }
1424
1425 EMData *result = new EMData();
1426 result->set_size(nx, ny, 1);
1427 result->to_zero();
1428 float *result_data = result->get_data();
1429
1430 EMData *this_copy = this;
1431 this_copy = copy();
1432
1433 for (int i = 0; i < nx; i++) {
1434 Transform t(Dict("type","2d","alpha",(float) M_PI * 2.0f * i / nx));
1435 this_copy->transform(t);
1436
1437 float *copy_data = this_copy->get_data();
1438
1439 for (int y = 0; y < nx; y++) {
1440 for (int x = 0; x < nx; x++) {
1441 if (Util::square(x - nx / 2) + Util::square(y - nx / 2) <= nx * nx / 4) {
1442 result_data[i + y * nx] += copy_data[x + y * nx];
1443 }
1444 }
1445 }
1446
1447 this_copy->update();
1448 }
1449
1450 result->update();
1451
1452 if( this_copy )
1453 {
1454 delete this_copy;
1455 this_copy = 0;
1456 }
1457
1458 EXITFUNC;
1459 return result;
1460}
1461
1462void EMData::zero_corner_circulant(const int radius)
1463{
1464 if ( nz > 1 && nz < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - nz is too small");
1465 if ( ny > 1 && ny < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - ny is too small");
1466 if ( nx > 1 && nx < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - nx is too small");
1467
1468 int it_z = radius;
1469 int it_y = radius;
1470 int it_x = radius;
1471
1472 if ( nz == 1 ) it_z = 0;
1473 if ( ny == 1 ) it_y = 0;
1474 if ( nx == 1 ) it_z = 0;
1475
1476 if ( nz == 1 && ny == 1 )
1477 {
1478 for ( int x = -it_x; x <= it_x; ++x )
1479 get_value_at_wrap(x) = 0;
1480
1481 }
1482 else if ( nz == 1 )
1483 {
1484 for ( int y = -it_y; y <= it_y; ++y)
1485 for ( int x = -it_x; x <= it_x; ++x )
1486 get_value_at_wrap(x,y) = 0;
1487 }
1488 else
1489 {
1490 for( int z = -it_z; z <= it_z; ++z )
1491 for ( int y = -it_y; y <= it_y; ++y)
1492 for ( int x = -it_x; x < it_x; ++x )
1493 get_value_at_wrap(x,y,z) = 0;
1494
1495 }
1496
1497}
1498
1499EMData *EMData::calc_ccf(EMData * with, fp_flag fpflag,bool center)
1500{
1501 ENTERFUNC;
1502
1503 if( with == 0 ) {
1504#ifdef EMAN2_USING_CUDA //CUDA
1505 if(EMData::usecuda == 1 && cudarwdata) {
1506 //cout << "calc ccf" << endl;
1507 EMData* ifft = 0;
1508 bool delifft = false;
1509 int offset = 0;
1510
1511 //do fft if not alreay done
1512 if(!is_complex()){
1513 ifft = do_fft_cuda();
1514 delifft = true;
1515 offset = 2 - nx%2;
1516 }else{
1517 ifft = this;
1518 }
1519 calc_conv_cuda(ifft->cudarwdata,ifft->cudarwdata,nx + offset, ny, nz); //this is the business end, results go in afft
1520
1521 EMData * conv = ifft->do_ift_cuda();
1522 if(delifft) delete ifft;
1523 conv->update();
1524
1525 return conv;
1526 }
1527#endif
1528 EXITFUNC;
1529 return convolution(this,this,fpflag, center);
1530 }
1531 else if ( with == this ){ // this if statement is not necessary, the correlation function tests to see if with == this
1532 EXITFUNC;
1533 return correlation(this, this, fpflag,center);
1534 }
1535 else {
1536
1537#ifdef EMAN2_USING_CUDA //CUDA
1538 // assume always get rw data (makes life a lot easier!!!
1539 // also assume that both images are the same size. When using CUDA we are only interested in speed, not flexibility!!
1540 // P.S. (I feel like I am pounding square pegs into a round holes with CUDA)
1541 if(EMData::usecuda == 1 && cudarwdata && with->cudarwdata) {
1542 //cout << "using CUDA for ccf" << endl;
1543 EMData* afft = 0;
1544 EMData* bfft = 0;
1545 bool delafft = false, delbfft = false;
1546 int offset = 0;
1547
1548 //do ffts if not alreay done
1549 if(!is_complex()){
1550 afft = do_fft_cuda();
1551 delafft = true;
1552 offset = 2 - nx%2;
1553 //cout << "Do cuda FFT A" << endl;
1554 }else{
1555 afft = this;
1556 }
1557 if(!with->is_complex()){
1558 bfft = with->do_fft_cuda();
1559 //cout << "Do cuda FFT B" << endl;
1560 delbfft = true;
1561 }else{
1562 bfft = with;
1563 }
1564
1565 calc_ccf_cuda(afft->cudarwdata,bfft->cudarwdata,nx + offset, ny, nz); //this is the business end, results go in afft
1566
1567 if(delbfft) delete bfft;
1568
1569 EMData * corr = afft->do_ift_cuda();
1570 if(delafft) delete afft;
1571 //cor->do_ift_inplace_cuda();//a bit faster, but I'll alos need to rearrnage the mem structure for it to work, BUT this is very SLOW.
1572 corr->update();
1573
1574 return corr;
1575 }
1576#endif
1577
1578 // If the argument EMData pointer is not the same size we automatically resize it
1579 bool undoresize = false;
1580 int wnx = with->get_xsize(); int wny = with->get_ysize(); int wnz = with->get_zsize();
1581 if (!(is_complex()^with->is_complex()) && (wnx != nx || wny != ny || wnz != nz) ) {
1582 Region r((wnx-nx)/2, (wny-ny)/2, (wnz-nz)/2,nx,ny,nz);
1583 with->clip_inplace(r);
1584 undoresize = true;
1585 }
1586
1587 EMData* cor = correlation(this, with, fpflag, center);
1588
1589 // If the argument EMData pointer was resized, it is returned to its original dimensions
1590 if ( undoresize ) {
1591 Region r((nx-wnx)/2, (ny-wny)/2,(nz-wnz)/2,wnx,wny,wnz);
1592 with->clip_inplace(r);
1593 }
1594
1595 EXITFUNC;
1596 return cor;
1597 }
1598}
1599
1600EMData *EMData::calc_ccf_masked(EMData *with,EMData *withsq,EMData *mask)
1601{
1602 if ((withsq==0 && mask!=0)||(withsq!=0 && mask==0))
1603 throw NullPointerException("calc_ccf_masked error, both or neither of withsq and mask must be specified");
1604
1605 bool needfree=0;
1606 if (withsq==0) {
1607 withsq=with->process("math.squared");
1608 mask=this->process("threshold.notzero");
1609 needfree=1;
1610 }
1611
1612 EMData *c1=this->calc_ccf(with);
1613 EMData *c2=mask->calc_ccf(withsq);
1614
1615
1616 c2->process_inplace("math.reciprocal",Dict("zero_to",0.0f));
1617 c2->process_inplace("math.sqrt");
1618 c1->mult(*c2);
1619 delete c2;
1620
1621 if (needfree) {
1622 delete withsq;
1623 delete mask;
1624 }
1625 return c1;
1626}
1627
1628EMData *EMData::calc_ccfx( EMData * const with, int y0, int y1, bool no_sum, bool flip,bool usez)
1629{
1630 ENTERFUNC;
1631
1632 if (!with) {
1633 LOGERR("NULL 'with' image. ");
1634 throw NullPointerException("NULL input image");
1635 }
1636
1637 if (!EMUtil::is_same_size(this, with)) {
1638 LOGERR("images not same size: (%d,%d,%d) != (%d,%d,%d)",
1639 nx, ny, nz,
1640 with->get_xsize(), with->get_ysize(), with->get_zsize());
1641 throw ImageFormatException("images not same size");
1642 }
1643 if (get_ndim() > 2) {
1644 LOGERR("2D images only");
1645 throw ImageDimensionException("2D images only");
1646 }
1647
1648 if (y1 <= y0) {
1649 y1 = ny;
1650 }
1651
1652 if (y0 >= y1) {
1653 y0 = 0;
1654 }
1655
1656 if (y0 < 0) {
1657 y0 = 0;
1658 }
1659
1660 if (y1 > ny) {
1661 y1 = ny;
1662 }
1663 if (is_complex_x() || with->is_complex_x() ) throw; // Woops don't support this anymore!
1664
1665// static int nx_fft = 0;
1666// static int ny_fft = 0;
1667// static EMData f1;
1668// static EMData f2;
1669// static EMData rslt;
1670
1671 int height = y1-y0;
1672 int width = (nx+2-(nx%2));
1673 int wpad = ((width+3)/4)*4; // This is for 128 bit alignment of rows to prevent SSE crashes
1674
1675 EMData f1(wpad,height);
1676 EMData f2(wpad,height);
1677 EMData rslt(nx,height);
1678
1679// if (wpad != nx_fft || height != ny_fft ) { // Seems meaningless, but due to static definitions above. f1,f2 are cached to prevent multiple reallocations
1680// f1.set_size(wpad,height);
1681// f2.set_size(wpad,height);
1682// rslt.set_size(nx,height);
1683// nx_fft = wpad;
1684// ny_fft = height;
1685// }
1686
1687#ifdef EMAN2_USING_CUDA
1688 // FIXME : Not tested with new wpad change
1689 if (EMData::usecuda == 1 && cudarwdata && with->cudarwdata) {
1690 //cout << "calc_ccfx CUDA" << endl;
1691 if(!f1.cudarwdata) f1.rw_alloc();
1692 if(!f2.cudarwdata) f2.rw_alloc();
1693 if(!rslt.cudarwdata) rslt.rw_alloc();
1694 cuda_dd_fft_real_to_complex_nd(cudarwdata, f1.cudarwdata, nx, 1, 1, height);
1695 cuda_dd_fft_real_to_complex_nd(with->cudarwdata, f2.cudarwdata, nx, 1, 1, height);
1696 calc_ccf_cuda(f1.cudarwdata, f2.cudarwdata, nx, ny, nz);
1697 cuda_dd_fft_complex_to_real_nd(f1.cudarwdata, rslt.cudarwdata, nx, 1, 1, height);
1698 if(no_sum){
1699 EMData* result = new EMData(rslt);
1700 return result;
1701 }
1702 EMData* cf = new EMData(0,0,nx,1,1); //cuda constructor
1703 cf->runcuda(emdata_column_sum(rslt.cudarwdata, nx, ny));
1704 cf->update();
1705
1706 EXITFUNC;
1707 return cf;
1708 }
1709#endif
1710
1711// printf("%d %d %d\n",(int)get_attr("nx"),(int)f2.get_attr("nx"),width);
1712
1713 float *d1 = get_data();
1714 float *d2 = with->get_data();
1715 float *f1d = f1.get_data();
1716 float *f2d = f2.get_data();
1717 for (int j = 0; j < height; j++) {
1718 EMfft::real_to_complex_1d(d1 + j * nx, f1d+j*wpad, nx);
1719 EMfft::real_to_complex_1d(d2 + j * nx, f2d+j*wpad, nx);
1720 }
1721
1722 if(flip == false) {
1723 for (int j = 0; j < height; j++) {
1724 float *f1a = f1d + j * wpad;
1725 float *f2a = f2d + j * wpad;
1726
1727 for (int i = 0; i < width / 2; i++) {
1728 float re1 = f1a[2*i];
1729 float re2 = f2a[2*i];
1730 float im1 = f1a[2*i+1];
1731 float im2 = f2a[2*i+1];
1732
1733 f1d[j*wpad+i*2] = re1 * re2 + im1 * im2;
1734 f1d[j*wpad+i*2+1] = im1 * re2 - re1 * im2;
1735 }
1736 }
1737 } else {
1738 for (int j = 0; j < height; j++) {
1739 float *f1a = f1d + j * wpad;
1740 float *f2a = f2d + j * wpad;
1741
1742 for (int i = 0; i < width / 2; i++) {
1743 float re1 = f1a[2*i];
1744 float re2 = f2a[2*i];
1745 float im1 = f1a[2*i+1];
1746 float im2 = f2a[2*i+1];
1747
1748 f1d[j*wpad+i*2] = re1 * re2 - im1 * im2;
1749 f1d[j*wpad+i*2+1] = im1 * re2 + re1 * im2;
1750 }
1751 }
1752 }
1753
1754 float* rd = rslt.get_data();
1755 for (int j = y0; j < y1; j++) {
1756 EMfft::complex_to_real_1d(f1d+j*wpad, rd+j*nx, nx);
1757 }
1758
1759 // This converts the CCF values to Z values (in terms of standard deviations above the mean), on a per-row basis
1760 // The theory is that this should better weight the radii that contribute most strongly to the orientation determination
1761 if (usez) {
1762 for (int y=0; y<height; y++) {
1763 float mn=0,sg=0;
1764 for (int x=0; x<nx; x++) {
1765 mn+=rd[x+y*nx];
1766 sg+=rd[x+y*nx]*rd[x+y*nx];
1767 }
1768 mn/=(float)nx; //mean
1769 sg=std::sqrt(sg/(float)nx-mn*mn); //sigma
1770 if (sg==0) sg=1.0;
1771
1772 for (int x=0; x<nx; x++) rd[x+y*nx]=(rd[x+y*nx]-mn)/sg;
1773 }
1774 }
1775
1776 if (no_sum) {
1777 rslt.update(); // This is important in terms of the copy - the returned object won't have the correct flags unless we do this
1778 EXITFUNC;
1779 return new EMData(rslt);
1780 } else {
1781 EMData *cf = new EMData(nx,1,1);
1782 cf->to_zero();
1783 float *c = cf->get_data();
1784 for (int j = 0; j < height; j++) {
1785 for(int i = 0; i < nx; ++i) {
1786 c[i] += rd[i+j*nx];
1787 }
1788 }
1789 cf->update();
1790 EXITFUNC;
1791 return cf;
1792 }
1793}
1794
1795EMData *EMData::make_rotational_footprint_cmc( bool unwrap) {
1796 ENTERFUNC;
1797 update_stat();
1798 // Note that rotational_footprint caching saves a large amount of time
1799 // but this is at the expense of memory. Note that a policy is hardcoded here,
1800 // that is that caching is only employed when premasked is false and unwrap
1801 // is true - this is probably going to be what is used in most scenarios
1802 // as advised by Steve Ludtke - In terms of performance this caching doubles the metric
1803 // generated by e2speedtest.
1804 if ( rot_fp != 0 && unwrap == true) {
1805 return new EMData(*rot_fp);
1806 }
1807
1808 static EMData obj_filt;
1809 static int updsize = 0; // for threadsafety
1810 EMData* filt = &obj_filt;
1811 filt->set_complex(true);
1812 int delfilt = 0;
1813
1814 // The filter object is nothing more than a cached high pass filter
1815 // Ultimately it is used an argument to the EMData::mult(EMData,prevent_complex_multiplication (bool))
1816 // function in calc_mutual_correlation. Note that in the function the prevent_complex_multiplication
1817 // set to true, which is used for speed reasons.
1818 if (filt->get_xsize() != nx+2-(nx%2) || filt->get_ysize() != ny ||
1819 filt->get_zsize() != nz ) {
1820 // In case of a threading conflict, this will ignore the cache for this specific computation
1821 if (updsize) {
1822 delfilt=1;
1823 printf("FAILSAFE!!!!\n");
1824 filt=new EMData(nx+2-(nx%2), ny, nz);
1825 filt->to_one();
1826 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
1827 }
1828 else {
1829 updsize=1;
1830 filt->set_size(nx+2-(nx%2), ny, nz);
1831 filt->to_one();
1832
1833 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
1834 }
1835 }
1836
1837 EMData *ccf = this->calc_mutual_correlation(this, true,filt);
1838
1839 ccf->sub(ccf->get_edge_mean());
1840 EMData *result = ccf->unwrap();
1841 delete ccf; ccf = 0;
1842 updsize=0;
1843 if (delfilt) delete filt;
1844
1845 EXITFUNC;
1846 if ( unwrap == true)
1847 {
1848 // this if statement reflects a strict policy of caching in only one scenario see comments at beginning of function block
1849
1850// Note that the if statement at the beginning of this function ensures that rot_fp is not zero, so there is no need
1851// to throw any exception
1852// if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
1853
1854// Here is where the caching occurs - the rot_fp takes ownsherhip of the pointer, and a deep copied EMData object is returned.
1855// The deep copy invokes a cost in terms of CPU cycles and memory, but prevents the need for complicated memory management (reference counting)
1856 rot_fp = result;
1857 return new EMData(*rot_fp);
1858 }
1859 else return result;
1860}
1861
1862EMData *EMData::make_rotational_footprint( bool unwrap) {
1863 ENTERFUNC;
1864 update_stat();
1865 // Note that rotational_footprint caching saves a large amount of time
1866 // but this is at the expense of memory. Note that a policy is hardcoded here,
1867 // that is that caching is only employed when premasked is false and unwrap
1868 // is true - this is probably going to be what is used in most scenarios
1869 // as advised by Steve Ludtke - In terms of performance this caching doubles the metric
1870 // generated by e2speedtest.
1871 if ( rot_fp != 0 && unwrap == true) {
1872 return new EMData(*rot_fp);
1873 }
1874
1875 EMData* ccf = this->calc_ccf(this,CIRCULANT,true);
1876// EMData* ccf = this->calc_ccf(this,PADDED,true); # this would probably be a bit better, but takes 4x longer :^/
1877 ccf->sub(ccf->get_edge_mean());
1878 //ccf->process_inplace("xform.phaseorigin.tocenter"); ccf did the centering
1879 EMData *result = ccf->unwrap();
1880 delete ccf; ccf = 0;
1881
1882 EXITFUNC;
1883 if ( unwrap == true)
1884 { // this if statement reflects a strict policy of caching in only one scenario see comments at beginning of function block
1885
1886// Note that the if statement at the beginning of this function ensures that rot_fp is not zero, so there is no need
1887// to throw any exception
1888// if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
1889
1890// Here is where the caching occurs - the rot_fp takes ownsherhip of the pointer, and a deep copied EMData object is returned.
1891// The deep copy invokes a cost in terms of CPU cycles and memory, but prevents the need for complicated memory management (reference counting)
1892 rot_fp = result;
1893 return new EMData(*rot_fp);
1894 }
1895 else return result;
1896}
1897
1898EMData *EMData::make_rotational_footprint_e1( bool unwrap)
1899{
1900 ENTERFUNC;
1901
1902 update_stat();
1903 // Note that rotational_footprint caching saves a large amount of time
1904 // but this is at the expense of memory. Note that a policy is hardcoded here,
1905 // that is that caching is only employed when premasked is false and unwrap
1906 // is true - this is probably going to be what is used in most scenarios
1907 // as advised by Steve Ludtke - In terms of performance this caching doubles the metric
1908 // generated by e2speedtest.
1909 if ( rot_fp != 0 && unwrap == true) {
1910 return new EMData(*rot_fp);
1911 }
1912
1913 static EMData obj_filt;
1914 EMData* filt = &obj_filt;
1915 filt->set_complex(true);
1916// Region filt_region;
1917
1918// if (nx & 1) {
1919// LOGERR("even image xsize only"); throw ImageFormatException("even image xsize only");
1920// }
1921
1922 int cs = (((nx * 7 / 4) & 0xfffff8) - nx) / 2; // this pads the image to 1 3/4 * size with result divis. by 8
1923
1924 static EMData big_clip;
1925 int big_x = nx+2*cs;
1926 int big_y = ny+2*cs;
1927 int big_z = 1;
1928 if ( nz != 1 ) {
1929 big_z = nz+2*cs;
1930 }
1931
1932
1933 if ( big_clip.get_xsize() != big_x || big_clip.get_ysize() != big_y || big_clip.get_zsize() != big_z ) {
1934 big_clip.set_size(big_x,big_y,big_z);
1935 }
1936 // It is important to set all newly established pixels around the boundaries to the mean
1937 // If this is not done then the associated rotational alignment routine breaks, in fact
1938 // everythin just goes foo.
1939
1940 big_clip.to_value(get_edge_mean());
1941
1942 if (nz != 1) {
1943 big_clip.insert_clip(this,IntPoint(cs,cs,cs));
1944 } else {
1945 big_clip.insert_clip(this,IntPoint(cs,cs,0));
1946 }
1947
1948 // The filter object is nothing more than a cached high pass filter
1949 // Ultimately it is used an argument to the EMData::mult(EMData,prevent_complex_multiplication (bool))
1950 // function in calc_mutual_correlation. Note that in the function the prevent_complex_multiplication
1951 // set to true, which is used for speed reasons.
1952 if (filt->get_xsize() != big_clip.get_xsize() +2-(big_clip.get_xsize()%2) || filt->get_ysize() != big_clip.get_ysize() ||
1953 filt->get_zsize() != big_clip.get_zsize()) {
1954 filt->set_size(big_clip.get_xsize() + 2-(big_clip.get_xsize()%2), big_clip.get_ysize(), big_clip.get_zsize());
1955 filt->to_one();
1956 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
1957#ifdef EMAN2_USING_CUDA
1958 /*
1959 if(EMData::usecuda == 1 && big_clip.cudarwdata)
1960 {
1961 filt->copy_to_cuda(); // since this occurs just once for many images, we don't pay much of a speed pentalty here, and we avoid the hassel of messing with sparx
1962 }
1963 */
1964#endif
1965 }
1966#ifdef EMAN2_USING_CUDA
1967 /*
1968 if(EMData::usecuda == 1 && big_clip.cudarwdata && !filt->cudarwdata)
1969 {
1970 filt->copy_to_cuda(); // since this occurs just once for many images, we don't pay much of a speed pentalty here, and we avoid the hassel of messing with sparx
1971 }
1972 */
1973#endif
1974
1975 EMData *mc = big_clip.calc_mutual_correlation(&big_clip, true,filt);
1976 mc->sub(mc->get_edge_mean());
1977
1978 static EMData sml_clip;
1979 int sml_x = nx * 3 / 2;
1980 int sml_y = ny * 3 / 2;
1981 int sml_z = 1;
1982 if ( nz != 1 ) {
1983 sml_z = nz * 3 / 2;
1984 }
1985
1986 if ( sml_clip.get_xsize() != sml_x || sml_clip.get_ysize() != sml_y || sml_clip.get_zsize() != sml_z ) {
1987 sml_clip.set_size(sml_x,sml_y,sml_z); }
1988 if (nz != 1) {
1989 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,-cs+nz/4));
1990 } else {
1991 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,0)); // same as padding change above
1992 }
1993
1994 delete mc; mc = 0;
1995 EMData * result = NULL;
1996
1997 if (nz == 1) {
1998 if (!unwrap) {
1999#ifdef EMAN2_USING_CUDA
2000 //if(EMData::usecuda == 1 && sml_clip.cudarwdata) throw UnexpectedBehaviorException("shap masking is not yet supported by CUDA");
2001#endif
2002 result = sml_clip.process("mask.sharp", Dict("outer_radius", -1, "value", 0));
2003
2004 }
2005 else {
2006 result = sml_clip.unwrap();
2007 }
2008 }
2009 else {
2010 // I am not sure why there is any consideration of non 2D images, but it was here
2011 // in the first port so I kept when I cleaned this function up (d.woolford)
2012// result = clipped_mc;
2013 result = new EMData(sml_clip);
2014 }
2015
2016#ifdef EMAN2_USING_CUDA
2017 //if (EMData::usecuda == 1) sml_clip.roneedsanupdate(); //If we didn't do this then unwrap would use data from the previous call of this function, happens b/c sml_clip is static
2018#endif
2019 EXITFUNC;
2020 if ( unwrap == true)
2021 { // this if statement reflects a strict policy of caching in only one scenario see comments at beginning of function block
2022
2023 // Note that the if statement at the beginning of this function ensures that rot_fp is not zero, so there is no need
2024 // to throw any exception
2025 if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
2026
2027 // Here is where the caching occurs - the rot_fp takes ownsherhip of the pointer, and a deep copied EMData object is returned.
2028 // The deep copy invokes a cost in terms of CPU cycles and memory, but prevents the need for complicated memory management (reference counting)
2029 rot_fp = result;
2030 return new EMData(*rot_fp);
2031 }
2032 else return result;
2033}
2034
2035EMData *EMData::make_footprint(int type)
2036{
2037// printf("Make fp %d\n",type);
2038 if (type==0) {
2039 EMData *un=make_rotational_footprint_e1(); // Use EMAN1's footprint strategy
2040 if (un->get_ysize() <= 6) {
2041 throw UnexpectedBehaviorException("In EMData::make_footprint. The rotational footprint is too small");
2042 }
2043 EMData *tmp=un->get_clip(Region(0,4,un->get_xsize(),un->get_ysize()-6)); // 4 and 6 are empirical
2044 EMData *cx=tmp->calc_ccfx(tmp,0,-1,1);
2045 EMData *fp=cx->get_clip(Region(0,0,cx->get_xsize()/2,cx->get_ysize()));
2046 delete un;
2047 delete tmp;
2048 delete cx;
2049 return fp;
2050 }
2051 else if (type==1 || type==2 ||type==5 || type==6) {
2052 int i,j,kx,ky,lx,ly;
2053
2054 EMData *fft=do_fft();
2055
2056 // map for x,y -> radius for speed
2057 int rmax=(get_xsize()+1)/2;
2058 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
2059 for (i=0; i<rmax; i++) {
2060 for (j=0; j<rmax; j++) {
2061 rmap[i+j*rmax]=hypot((float)i,(float)j);
2062// printf("%d\t%d\t%f\n",i,j,rmap[i+j*rmax]);
2063 }
2064 }
2065
2066 EMData *fp=new EMData(rmax*2+2,rmax*2,1);
2067 fp->set_complex(1);
2068 fp->to_zero();
2069
2070 // Two vectors in to complex space (kx,ky) and (lx,ly)
2071 // We are computing the bispectrum, f(k).f(l).f*(k+l)
2072 // but integrating out two dimensions, leaving |k|,|l|
2073 for (kx=-rmax+1; kx<rmax; kx++) {
2074 for (ky=-rmax+1; ky<rmax; ky++) {
2075 for (lx=-rmax+1; lx<rmax; lx++) {
2076 for (ly=-rmax+1; ly<rmax; ly++) {
2077 int ax=kx+lx;
2078 int ay=ky+ly;
2079 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
2080 int r1=(int)floor(.5+rmap[abs(kx)+rmax*abs(ky)]);
2081 int r2=(int)floor(.5+rmap[abs(lx)+rmax*abs(ly)]);
2082// if (r1>500 ||r2>500) printf("%d\t%d\t%d\t%d\t%d\t%d\n",kx,ky,lx,ly,r1,r2);
2083// float r3=rmap[ax+rmax*ay];
2084 if (r1+r2>=rmax) continue;
2085
2086 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
2087 fp->set_value_at(r1*2,r2,p.real()+fp->get_value_at(r1*2,r2)); // We keep only the real component in anticipation of zero phase sum
2088// fp->set_value_at(r1*2,rmax*2-r2-1, fp->get_value_at(r1*2,r2)); // We keep only the real component in anticipation of zero phase sum
2089// fp->set_value_at(r1*2+1,r2,p.real()+fp->get_value_at(r1*2+1,r2)); // We keep only the real component in anticipation of zero phase sum
2090 fp->set_value_at(r1*2+1,r2,fp->get_value_at(r1*2+1,r2)+1); // a normalization counter
2091 }
2092 }
2093 }
2094 }
2095
2096 // Normalizes the pixels based on the accumulated counts then sets the imaginary components back to zero
2097 if (type==5 || type==6) {
2098 for (i=0; i<rmax*2; i+=2) {
2099 for (j=0; j<rmax; j++) {
2100 float norm=fp->get_value_at(i+1,j);
2101 fp->set_value_at(i,rmax*2-j-1,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
2102 fp->set_value_at(i,j,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
2103 fp->set_value_at(i+1,j,0.0);
2104 }
2105 }
2106 }
2107 else {
2108 for (i=0; i<rmax*2; i+=2) {
2109 for (j=0; j<rmax; j++) {
2110 float norm=fp->get_value_at(i+1,j);
2111 fp->set_value_at(i,rmax*2-j-1,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
2112 fp->set_value_at(i,j,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
2113 fp->set_value_at(i+1,j,0.0);
2114 }
2115 }
2116 }
2117
2118 free(rmap);
2119 if (type==2||type==6) {
2120 EMData *f2=fp->do_ift();
2121 if (f2->get_value_at(0,0)<0) f2->mult(-1.0f);
2122 f2->process_inplace("xform.phaseorigin.tocorner");
2123 delete fp;
2124 return f2;
2125 }
2126 return fp;
2127 }
2128 else if (type==3 || type==4) {
2129 int h,i,j,kx,ky,lx,ly;
2130
2131 EMData *fft=do_fft();
2132
2133 // map for x,y -> radius for speed
2134 int rmax=(get_xsize()+1)/2;
2135 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
2136 for (i=0; i<rmax; i++) {
2137 for (j=0; j<rmax; j++) {
2138 rmap[i+j*rmax]=hypot((float)i,(float)j);
2139// printf("%d\t%d\t%f\n",i,j,rmap[i+j*rmax]);
2140 }
2141 }
2142
2143 EMData *fp=new EMData(rmax*2+2,rmax*2,16);
2144
2145 fp->set_complex(1);
2146 fp->to_zero();
2147
2148 // Two vectors in to complex space (kx,ky) and (lx,ly)
2149 // We are computing the bispectrum, f(k).f(l).f*(k+l)
2150 // but integrating out two dimensions, leaving |k|,|l|
2151 for (kx=-rmax+1; kx<rmax; kx++) {
2152 for (ky=-rmax+1; ky<rmax; ky++) {
2153 for (lx=-rmax+1; lx<rmax; lx++) {
2154 for (ly=-rmax+1; ly<rmax; ly++) {
2155 int ax=kx+lx;
2156 int ay=ky+ly;
2157 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
2158 float rr1=rmap[abs(kx)+rmax*abs(ky)];
2159 float rr2=rmap[abs(lx)+rmax*abs(ly)];
2160 int r1=(int)floor(.5+rr1);
2161 int r2=(int)floor(.5+rr2);
2162// if (r1>500 ||r2>500) printf("%d\t%d\t%d\t%d\t%d\t%d\n",kx,ky,lx,ly,r1,r2);
2163// float r3=rmap[ax+rmax*ay];
2164 if (r1+r2>=rmax || rr1==0 ||rr2==0) continue;
2165
2166 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
2167 int dot=(int)floor((kx*lx+ky*ly)/(rr1*rr2)*7.5); // projection of k on l 0-31
2168 if (dot<0) dot=16+dot;
2169// int dot=(int)floor((kx*lx+ky*ly)/(rr1*rr2)*7.5+8.0); // projection of k on l 0-15
2170 fp->set_value_at(r1*2,r2,dot,p.real()+fp->get_value_at(r1*2,r2,dot)); // We keep only the real component in anticipation of zero phase sum
2171// fp->set_value_at(r1*2,rmax*2-r2-1, fp->get_value_at(r1*2,r2)); // We keep only the real component in anticipation of zero phase sum
2172// fp->set_value_at(r1*2+1,r2,p.real()+fp->get_value_at(r1*2+1,r2)); // We keep only the real component in anticipation of zero phase sum
2173 fp->set_value_at(r1*2+1,r2,dot,fp->get_value_at(r1*2+1,r2,dot)+1); // a normalization counter
2174 }
2175 }
2176 }
2177 }
2178
2179 // Normalizes the pixels based on the accumulated counts then sets the imaginary components back to zero
2180 for (i=0; i<rmax*2; i+=2) {
2181 for (j=0; j<rmax; j++) {
2182 for (h=0; h<16; h++) {
2183 float norm=fp->get_value_at(i+1,j,h);
2184// fp->set_value_at(i,rmax*2-j-1,h,cbrt(fp->get_value_at(i,j,h)/(norm==0?1.0:norm)));
2185// fp->set_value_at(i,j,h,cbrt(fp->get_value_at(i,j,h)/(norm==0?1.0:norm)));
2186 fp->set_value_at(i,rmax*2-j-1,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
2187 fp->set_value_at(i,j,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
2188 // fp->set_value_at(i,rmax*2-j-1,fp->get_value_at(i,j)/(norm==0?1.0:norm));
2189 // fp->set_value_at(i,j,fp->get_value_at(i,j)/(norm==0?1.0:norm));
2190 fp->set_value_at(i+1,j,h,0.0);
2191 }
2192 }
2193 }
2194
2195 free(rmap);
2196 if (type==4) {
2197 EMData *f2=fp->do_ift();
2198 if (f2->get_value_at(0,0,0)<0) f2->mult(-1.0f);
2199 f2->process_inplace("xform.phaseorigin.tocorner");
2200 delete fp;
2201 return f2;
2202 }
2203 return fp;
2204 }
2205 throw UnexpectedBehaviorException("There is not implementation for the parameters you specified");
2206}
2207
2208
2209EMData *EMData::calc_mutual_correlation(EMData * with, bool tocenter, EMData * filter)
2210{
2211 ENTERFUNC;
2212
2213 if (with && !EMUtil::is_same_size(this, with)) {
2214 LOGERR("images not same size");
2215 throw ImageFormatException( "images not same size");
2216 }
2217
2218#ifdef EMAN2_USING_CUDA
2219 if(EMData::usecuda == 1 && cudarwdata && with->cudarwdata)
2220 {
2221
2222 EMData* this_fft = do_fft_cuda();
2223
2224 EMData *cf = 0;
2225 if (with && with != this) {
2226 cf = with->do_fft_cuda();
2227 }else{
2228 cf = this_fft->copy();
2229 }
2230
2231 if (filter) {
2232 if (!EMUtil::is_same_size(filter, cf)) {
2233 LOGERR("improperly sized filter");
2234 throw ImageFormatException("improperly sized filter");
2235 }
2236 mult_complex_efficient_cuda(cf->cudarwdata, filter->cudarwdata, cf->get_xsize(), cf->get_ysize(), cf->get_zsize(), 1);
2237 mult_complex_efficient_cuda(this_fft->cudarwdata, filter->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize(), 1);
2238 }
2239
2240 mcf_cuda(this_fft->cudarwdata, cf->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize());
2241
2242 EMData *f2 = cf->do_ift_cuda();
2243
2244 if (tocenter) {
2245 f2->process_inplace("xform.phaseorigin.tocenter");
2246 }
2247
2248 if( cf )
2249 {
2250 delete cf;
2251 cf = 0;
2252 }
2253
2254 if( this_fft )
2255 {
2256 delete this_fft;
2257 this_fft = 0;
2258 }
2259
2260 f2->set_attr("label", "MCF");
2261 f2->set_path("/tmp/eman.mcf");
2262 f2->update();
2263
2264 EXITFUNC;
2265 return f2;
2266 }
2267#endif
2268
2269 EMData *this_fft = 0;
2270 this_fft = do_fft();
2271
2272 if (!this_fft) {
2273
2274 LOGERR("FFT returns NULL image");
2275 throw NullPointerException("FFT returns NULL image");
2276 }
2277
2278 this_fft->ap2ri(); //this is not needed!
2279 EMData *cf = 0;
2280
2281 if (with && with != this) {
2282 cf = with->do_fft();
2283 if (!cf) {
2284 LOGERR("FFT returns NULL image");
2285 throw NullPointerException("FFT returns NULL image");
2286 }
2287 cf->ap2ri(); //nor is this!
2288 }
2289 else {
2290 cf = this_fft->copy();
2291 }
2292
2293 if (filter) {
2294 if (!EMUtil::is_same_size(filter, cf)) {
2295 LOGERR("improperly sized filter");
2296 throw ImageFormatException("improperly sized filter");
2297 }
2298
2299 cf->mult_complex_efficient(*filter,true); //insanely this is required....
2300 this_fft->mult(*filter,true);
2301 //cf->mult_complex_efficient(*filter,7); // takes advantage of the fact that the filter is 1 everywhere but near the origin
2302 //this_fft->mult_complex_efficient(*filter,7);
2303 /*cf->mult_complex_efficient(*filter,5);
2304 this_fft->mult_complex_efficient(*filter,5);*/
2305 }
2306
2307 float *rdata1 = this_fft->get_data();
2308 float *rdata2 = cf->get_data();
2309 size_t this_fft_size = (size_t)this_fft->get_xsize() * this_fft->get_ysize() * this_fft->get_zsize();
2310
2311 if (with == this) {
2312 for (size_t i = 0; i < this_fft_size; i += 2) {
2313 rdata2[i] = std::sqrt(rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
2314 rdata2[i + 1] = 0;
2315 }
2316
2317 this_fft->update();
2318 cf->update();
2319 }
2320 else {
2321 for (size_t i = 0; i < this_fft_size; i += 2) {
2322 rdata2[i] = (rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
2323 rdata2[i + 1] = (rdata1[i + 1] * rdata2[i] - rdata1[i] * rdata2[i + 1]);
2324 }
2325
2326 //This seems like a bug, but it probably is never used....
2327 for (size_t i = 0; i < this_fft_size; i += 2) {
2328 float t = Util::square(rdata2[i]) + Util::square(rdata2[i + 1]);
2329 if (t != 0) {
2330 t = pow(t, 0.25f);
2331 rdata2[i] /= t;
2332 rdata2[i + 1] /= t;
2333 }
2334 }
2335 this_fft->update();
2336 cf->update();
2337 }
2338
2339 EMData *f2 = cf->do_ift();
2340
2341 if (tocenter) {
2342 f2->process_inplace("xform.phaseorigin.tocenter");
2343 }
2344
2345 if( cf )
2346 {
2347 delete cf;
2348 cf = 0;
2349 }
2350
2351 if( this_fft )
2352 {
2353 delete this_fft;
2354 this_fft = 0;
2355 }
2356
2357 f2->set_attr("label", "MCF");
2358 f2->set_path("/tmp/eman.mcf");
2359
2360 EXITFUNC;
2361 return f2;
2362}
2363
2364
2365vector < float > EMData::calc_hist(int hist_size, float histmin, float histmax,const float& brt, const float& cont)
2366{
2367 ENTERFUNC;
2368
2369// static size_t prime[] = { 1, 3, 7, 11, 17, 23, 37, 59, 127, 253, 511 };
2370
2371 size_t size = (size_t)nx * ny * nz;
2372 if (histmin == histmax) {
2373 histmin = get_attr("minimum");
2374 histmax = get_attr("maximum");
2375 }
2376 int n=hist_size;
2377 float w = (float)(n-1) / (histmax - histmin);
2378 float * data = get_data();
2379
2380 vector <float> hist(hist_size, 0.0);
2381
2382 if (cont != 1.0f || brt != 0) {
2383 for(size_t i=0; i<size; i++) {
2384 float val = cont*(data[i]+brt);
2385 int j = Util::round((val - histmin) * w);
2386
2387 // Outliers now go in the edge bins
2388 if (j>=n) j=n-1;
2389 if (j<0) j=0;
2390 hist[j]+=1;
2391 }
2392 }
2393 else {
2394 for(size_t i=0; i<size; i++) {
2395 float val = data[i];
2396 int j = Util::round((val - histmin) * w);
2397 if (j>=n) j=n-1;
2398 if (j<0) j=0;
2399 hist[j]+=1;
2400 }
2401 }
2402
2403// 20 years ago, it was expensive to compute histograms of large images
2404// so there was a complicated strategy to approximate them
2405// int p0 = 0;
2406// int p1 = 0;
2407// size_t size = (size_t)nx * ny * nz;
2408// if (size < 300000) {
2409// p0 = 0;
2410// p1 = 0;
2411// }
2412// else if (size < 2000000) {
2413// p0 = 2;
2414// p1 = 3;
2415// }
2416// else if (size < 8000000) {
2417// p0 = 4;
2418// p1 = 6;
2419// }
2420// else {
2421// p0 = 7;
2422// p1 = 9;
2423// }
2424//
2425// if (is_complex() && p0 > 0) {
2426// p0++;
2427// p1++;
2428// }
2429
2430// size_t di = 0;
2431// // float norm = 0;
2432// size_t n = hist.size();
2433//
2434// float * data = get_data();
2435// for (int k = p0; k <= p1; ++k) {
2436// if (is_complex()) {
2437// di = prime[k] * 2;
2438// }
2439// else {
2440// di = prime[k];
2441// }
2442//
2443// // norm += (float)size / (float) di;
2444// float w = (float)n / (histmax - histmin);
2445//
2446// for(size_t i=0; i<=size-di; i += di) {
2447// float val;
2448// if (cont != 1.0f || brt != 0)val = cont*(data[i]+brt);
2449// else val = data[i];
2450// int j = Util::round((val - histmin) * w);
2451// if (j >= 0 && j < (int) n) {
2452// hist[j] += 1;
2453// }
2454// }
2455// }
2456// /*
2457// for (size_t i = 0; i < hist.size(); ++i) {
2458// if (norm != 0) {
2459// hist[i] = hist[i] / norm;
2460// }
2461// }
2462// */
2463 return hist;
2464
2465 EXITFUNC;
2466}
2467
2468
2469
2470
2471
2472vector<float> EMData::calc_az_dist(int n, float a0, float da, float rmin, float rmax)
2473{
2474 ENTERFUNC;
2475
2476 a0=a0*M_PI/180.0f;
2477 da=da*M_PI/180.0f;
2478
2479 if (get_ndim() > 2) {
2480 throw ImageDimensionException("no 3D image");
2481 }
2482
2483 float *yc = new float[n];
2484
2485 vector<float> vd(n);
2486 for (int i = 0; i < n; i++) {
2487 vd[i]=0;
2488 yc[i] = 0.00001f;
2489 }
2490
2491 int isri=is_ri();
2492
2493 float * data = get_data();
2494 if (is_complex()) {
2495 int c = 0;
2496 for (int y = 0; y < ny; y++) {
2497 for (int x = 0; x < nx; x += 2, c += 2) {
2498 int x1 = x / 2;
2499 int y1 = y<ny/2?y:y-ny;
2500 float r = (float)Util::hypot_fast(x1,y1);
2501
2502 if (r >= rmin && r <= rmax) {
2503 float a = 0;
2504
2505 if (y != ny / 2 || x != 0) {
2506 a = (atan2((float)y1, (float)x1) - a0) / da;
2507 }
2508
2509 int i = (int)(floor(a));
2510 a -= i;
2511
2512 float h=0;
2513 if (isri) {
2514 h = (float)hypot(data[c], data[c + 1]);
2515 }
2516 else h = data[c];
2517 if (i == 0) {
2518 vd[0] += h * (1.0f - a);
2519 yc[0] += (1.0f - a);
2520 }
2521 else if (i == n - 1) {
2522 vd[n - 1] += h * a;
2523 yc[n - 1] += a;
2524 }
2525 else if (i > 0 && i < (n - 1)) {
2526 vd[i] += h * (1.0f - a);
2527 yc[i] += (1.0f - a);
2528 vd[i + 1] += h * a;
2529 yc[i + 1] += a;
2530 }
2531 }
2532 }
2533 }
2534 }
2535 else {
2536 int c = 0;
2537 float half_nx = (nx - 1) / 2.0f;
2538 float half_ny = (ny - 1) / 2.0f;
2539
2540 for (int y = 0; y < ny; y++) {
2541 for (int x = 0; x < nx; x++, c++) {
2542 float y1 = y - half_ny;
2543 float x1 = x - half_nx;
2544 float r = (float)hypot(x1, y1);
2545
2546 if (r >= rmin && r <= rmax) {
2547 float a = 0;
2548 if (x1 != 0 || y1 != 0) {
2549 a = atan2(y1, x1);
2550 if (a < 0) {
2551 a += M_PI * 2;
2552 }
2553 }
2554
2555 a = (a - a0) / da;
2556 int i = (int)(floor(a));
2557 a -= i;
2558
2559 if (i == 0) {
2560 vd[0] += data[c] * (1.0f - a);
2561 yc[0] += (1.0f - a);
2562 }
2563 else if (i == n - 1) {
2564 vd[n - 1] += data[c] * a;
2565 yc[n - 1] += (a);
2566 }
2567 else if (i > 0 && i < (n - 1)) {
2568 vd[i] += data[c] * (1.0f - a);
2569 yc[i] += (1.0f - a);
2570 vd[i + 1] += data[c] * a;
2571 yc[i + 1] += a;
2572 }
2573 }
2574 }
2575 }
2576 }
2577
2578
2579 for (int i = 0; i < n; i++) {
2580 if (vd[i]<0||yc[i]<0) printf("%d vd %f yc %f\n",i,vd[i],yc[i]);
2581 vd[i] /= yc[i];
2582 }
2583
2584 if( yc )
2585 {
2586 delete[]yc;
2587 yc = 0;
2588 }
2589
2590 return vd;
2591
2592 EXITFUNC;
2593}
2594
2595
2596EMData *EMData::unwrap(int r1, int r2, int xs, int dx, int dy, bool do360, bool weight_radial) const
2597{
2598 ENTERFUNC;
2599
2600 if (get_ndim() != 2) {
2601 throw ImageDimensionException("2D image only");
2602 }
2603
2604 int p = 1;
2605 if (do360) {
2606 p = 2;
2607 }
2608
2609 if (xs < 1) {
2610 xs = (int) Util::fast_floor(p * M_PI * ny / 3.0);
2611 xs-=xs%4; // 128 bit alignment, though best_fft_size may override
2612 xs = Util::calc_best_fft_size(xs);
2613 if (xs<=8) xs=16;
2614 }
2615
2616 if (r1 < 0) {
2617 r1 = 4;
2618 }
2619
2620 int rr = ny / 2 - 2 - (int) Util::fast_floor((float)(hypot(dx, dy)));
2621 rr-=rr%2;
2622 if (r2 <= r1 || r2 > rr) {
2623 r2 = rr;
2624 }
2625
2626 if ( (r2-r1) < 0 ) throw UnexpectedBehaviorException("The combination of function the arguments and the image dimensions causes unexpected behavior internally. Use a larger image, or a smaller value of r1, or a combination of both");
2627
2628#ifdef EMAN2_USING_CUDA
2629 if (EMData::usecuda == 1 && isrodataongpu()){
2630 //cout << " CUDA unwrap" << endl;
2631 EMData* result = new EMData(0,0,xs,(r2-r1),1);
2632 if(!result->rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");
2633 bindcudaarrayA(true);
2634 emdata_unwrap(result->cudarwdata, r1, r2, xs, p, dx, dy, weight_radial, nx, ny);
2635 unbindcudaarryA();
2636 result->update();
2637 return result;
2638 }
2639#endif
2640
2641 EMData *ret = new EMData();
2642 ret->set_size(xs, r2 - r1, 1);
2643 const float *const d = get_const_data();
2644 float *dd = ret->get_data();
2645 float pfac = (float)p/(float)xs;
2646 int nxon2 = nx/2;
2647 int nyon2 = ny/2;
2648 for (int x = 0; x < xs; x++) {
2649 float ang = x * M_PI * pfac;
2650 float si = sin(ang);
2651 float co = cos(ang);
2652
2653 for (int y = 0; y < r2 - r1; y++) {
2654 float ypr1 = (float)y + r1;
2655 float xx = ypr1 * co + nxon2 + dx;
2656 float yy = ypr1 * si + nyon2 + dy;
2657// float t = xx - Util::fast_floor(xx);
2658// float u = yy - Util::fast_floor(yy);
2659 float t = xx - (int)xx;
2660 float u = yy - (int)yy;
2661// int k = (int) Util::fast_floor(xx) + (int) (Util::fast_floor(yy)) * nx;
2662 int k = (int) xx + ((int) yy) * nx;
2663 float val = Util::bilinear_interpolate(d[k], d[k + 1], d[k + nx], d[k + nx+1], t,u);
2664 if (weight_radial) val *= ypr1;
2665 dd[x + y * xs] = val;
2666 }
2667
2668 }
2669 ret->update();
2670
2671 EXITFUNC;
2672 return ret;
2673}
2674
2675// NOTE : x axis is from 0 to 0.5 (Nyquist), and thus properly handles non-square images
2676// complex only
2677void EMData::apply_radial_func(float x0, float step, vector < float >array, bool interp)
2678{
2679 ENTERFUNC;
2680
2681 if (!is_complex()) throw ImageFormatException("apply_radial_func requires a complex image");
2682
2683 int n = (int)(array.size());
2684
2685 if (n*step>2.0) printf("Warning, apply_radial_func takes x0 and step with respect to Nyquist (0.5)\n");
2686
2687// printf("%f %f %f\n",array[0],array[25],array[50]);
2688
2689 ap2ri();
2690
2691 size_t ndims = get_ndim();
2692 float * data = get_data();
2693 if (ndims == 2) {
2694 int k = 0;
2695 for (int j = 0; j < ny; j++) {
2696 for (int i = 0; i < nx; i += 2, k += 2) {
2697 float r;
2698 if (j<ny/2) r = (float)hypot(i/(float)(nx*2), j/(float)ny);
2699 else r = (float)hypot(i/(float)(nx*2), (ny-j)/(float)ny);
2700 r = (r - x0) / step;
2701
2702 int l = 0;
2703 if (interp) {
2704 l = (int) floor(r);
2705 }
2706 else {
2707 l = (int) floor(r + 1);
2708 }
2709
2710
2711 float f = 0;
2712 if (l >= n - 2) {
2713 f = array[n - 1];
2714 }
2715 else {
2716 if (interp) {
2717 r -= l;
2718 f = (array[l] * (1.0f - r) + array[l + 1] * r);
2719 }
2720 else {
2721 f = array[l];
2722 }
2723 }
2724
2725 data[k] *= f;
2726 data[k + 1] *= f;
2727 }
2728 }
2729 }
2730 else if (ndims == 3) {
2731 int k = 0;
2732 for (int m = 0; m < nz; m++) {
2733 float mnz;
2734 if (m<nz/2) mnz=m*m/(float)(nz*nz);
2735 else { mnz=(nz-m)/(float)nz; mnz*=mnz; }
2736
2737 for (int j = 0; j < ny; j++) {
2738 float jny;
2739 if (j<ny/2) jny= j*j/(float)(ny*ny);
2740 else { jny=(ny-j)/(float)ny; jny*=jny; }
2741
2742 for (int i = 0; i < nx; i += 2, k += 2) {
2743 float r = std::sqrt((i * i / (nx*nx*4.0f)) + jny + mnz);
2744 r = (r - x0) / step;
2745
2746 int l = 0;
2747 if (interp) {
2748 l = (int) floor(r);
2749 }
2750 else {
2751 l = (int) floor(r + 1);
2752 }
2753
2754 float f = 0;
2755 if (l >= n - 2) {
2756 f = array[n - 1];
2757 }
2758 else {
2759 if (interp) {
2760 r -= l;
2761 f = (array[l] * (1.0f - r) + array[l + 1] * r);
2762 }
2763 else {
2764 f = array[l];
2765 }
2766 }
2767
2768//if (k%5000==0) printf("%d %d %d %f\n",i,j,m,f);
2769 data[k] *= f;
2770 data[k + 1] *= f;
2771 }
2772 }
2773 }
2774
2775 }
2776
2777 update();
2778 EXITFUNC;
2779}
2780
2781vector<float> EMData::calc_radial_dist(int n, float x0, float dx, int inten)
2782{
2783 ENTERFUNC;
2784
2785 vector<double>ret(n);
2786 vector<double>norm(n);
2787 vector<double>count(n);
2788
2789 int x,y,z,i;
2790 int step=is_complex()?2:1;
2791 int isinten=get_attr_default("is_intensity",0);
2792 int isri=is_ri();
2793
2794 if (isinten&&!inten) { throw InvalidParameterException("Must set inten for calc_radial_dist with intensity image"); }
2795
2796 switch (inten){
2797 case 0:
2798 case 1:
2799 case 4:
2800 case 5:
2801 for (i=0; i<n; i++) ret[i]=norm[i]=count[i]=0.0;
2802 break;
2803 case 2:
2804 for (i=0; i<n; i++) ret[i]=1.0e27;
2805 break;
2806 case 3:
2807 for (i=0; i<n; i++) ret[i]=-1.0e27;
2808 break;
2809 }
2810
2811 float * data = get_data();
2812
2813 // We do 2D separately to avoid the hypot3 call
2814 if (nz==1) {
2815 for (y=i=0; y<ny; y++) {
2816 for (x=0; x<nx; x+=step,i+=step) {
2817 float r,v;
2818 int f;
2819 if (step==2) { //complex
2820 if (x==0 && y>ny/2) continue;
2821 r=(float)(Util::hypot_fast(x/2,y<ny/2?y:ny-y)); // origin at 0,0; periodic
2822 r=(r-x0)/dx;
2823 f=int(r); // safe truncation, so floor isn't needed
2824 if (f<0 || f>=n) continue;
2825 switch (inten) {
2826 case 0:
2827 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2828 else v=data[i]; // amp/phase, just get amp
2829 break;
2830 case 1:
2831 if (isinten) v=data[i];
2832 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
2833 else v=data[i]*data[i];
2834 break;
2835 case 2:
2836 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2837 else v=data[i];
2838 if (v<ret[f]) ret[f]=v;
2839 break;
2840 case 3:
2841 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2842 else v=data[i];
2843 if (v>ret[f]) ret[f]=v;
2844 break;
2845 case 4:
2846 if (isinten) v=data[i];
2847 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
2848 else v=data[i]*data[i];
2849 ret[f]+=std::sqrt(v);
2850 norm[f]+=v;
2851 count[f]+=1.0;
2852 break;
2853 }
2854 }
2855 else {
2856 if (inten==5) r=(float)(Util::hypot_fast(x,y-ny/2));
2857 else r=(float)(Util::hypot_fast(x-nx/2,y-ny/2));
2858 r=(r-x0)/dx;
2859 f=int(r); // safe truncation, so floor isn't needed
2860 if (f<0 || f>=n) continue;
2861 switch (inten) {
2862 case 0:
2863 case 5:
2864 v=data[i];
2865 break;
2866 case 1:
2867 v=data[i]*data[i];
2868 break;
2869 case 2:
2870 if (data[i]<ret[f]) ret[f]=data[i];
2871 break;
2872 case 3:
2873 if (data[i]>ret[f]) ret[f]=data[i];
2874 break;
2875 case 4:
2876 ret[f]+=data[i];
2877 norm[f]+=data[i]*data[i];
2878 count[f]+=1.0;
2879 break;
2880 }
2881 }
2882
2883 if (inten<2||inten==5) {
2884 r-=float(f); // r is now the fractional spacing between bins
2885 // printf("%d\t%d\t%d\t%1.3f\t%d\t%1.3f\t%1.4g\n",x,y,f,r,step,Util::hypot_fast(x/2,y<ny/2?y:ny-y),v);
2886 ret[f]+=v*(1.0f-r);
2887 norm[f]+=(1.0f-r);
2888 if (f<n-1) {
2889 ret[f+1]+=v*r;
2890 norm[f+1]+=r;
2891 }
2892 }
2893 }
2894 }
2895 }
2896 else {
2897// FILE *out = fopen("x.txt","w");
2898 size_t i; //3D file may have >2G size
2899 for (z=i=0; z<nz; ++z) {
2900 for (y=0; y<ny; ++y) {
2901 for (x=0; x<nx; x+=step,i+=step) {
2902 float r,v;
2903 int f;
2904 if (step==2) { //complex
2905 if (x==0 && z>nz/2) continue;
2906 if (x==0 && z==nz/2 && y>ny/2) continue;
2907 r=Util::hypot3(x/2,y<ny/2?y:ny-y,z<nz/2?z:nz-z); // origin at 0,0; periodic
2908 r=(r-x0)/dx;
2909 f=int(r); // safe truncation, so floor isn't needed
2910 if (f<0 || f>=n) continue;
2911 switch(inten) {
2912 case 0:
2913 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2914 else v=data[i]; // amp/phase, just get amp
2915 break;
2916 case 1:
2917 if (isinten) v=data[i];
2918 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
2919 else v=data[i]*data[i];
2920 break;
2921 case 2:
2922 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2923 else v=data[i]; // amp/phase, just get amp
2924 if (v<ret[f]) ret[f]=v;
2925 break;
2926 case 3:
2927 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2928 else v=data[i]; // amp/phase, just get amp
2929 if (v>ret[f]) ret[f]=v;
2930 break;
2931 case 4:
2932 if (isinten) v=data[i];
2933 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
2934 else v=data[i]*data[i];
2935 ret[f]+=std::sqrt(v);
2936 norm[f]+=v;
2937 count[f]+=1.0;
2938// if (f==100) fprintf(out,"%1.4g\n",std::sqrt(v));
2939 break;
2940 }
2941 }
2942 else {
2943 if (inten==5) r=Util::hypot3(x,y-ny/2,z-nz/2);
2944 else r=Util::hypot3(x-nx/2,y-ny/2,z-nz/2);
2945 r=(r-x0)/dx;
2946 f=int(r); // safe truncation, so floor isn't needed
2947 if (f<0 || f>=n) continue;
2948 switch(inten) {
2949 case 0:
2950 case 5:
2951 v=data[i];
2952 break;
2953 case 1:
2954 v=data[i]*data[i];
2955 break;
2956 case 2:
2957 if (data[i]<ret[f]) ret[f]=data[i];
2958 break;
2959 case 3:
2960 if (data[i]>ret[f]) ret[f]=data[i];
2961 break;
2962 case 4:
2963 ret[f]+=data[i];
2964 norm[f]+=data[i]*data[i];
2965 count[f]+=1.0;
2966 break;
2967 }
2968 }
2969
2970 if (inten<2||inten==5) {
2971// ret[f]+=v;
2972// norm[f]+=1.0;
2973
2974 r-=float(f); // r is now the fractional spacing between bins
2975 ret[f]+=v*(1.0f-r);
2976 norm[f]+=(1.0f-r);
2977 if (f<n-1) {
2978 ret[f+1]+=v*r;
2979 norm[f+1]+=r;
2980 }
2981 }
2982 }
2983 }
2984 }
2985// fclose(out);
2986 }
2987
2988 if (inten<2||inten==5) {
2989 for (i=0; i<n; i++) ret[i]/=(norm[i]==0?1.0f:norm[i]); // Normalize
2990 }
2991 else if (inten==4) {
2992 for (i=0; i<n; i++) {
2993 ret[i]/=count[i]; // becomes mean
2994 norm[i]/=count[i]; // avg amp^2
2995 ret[i]*=ret[i];
2996 if (norm[i]<=ret[i]) ret[i]=0.0;
2997 else ret[i]=std::sqrt(norm[i]-ret[i]); // sigma
2998 }
2999 }
3000
3001 EXITFUNC;
3002
3003 return vector<float>(ret.begin(),ret.end());
3004}
3005
3006vector<float> EMData::calc_radial_dist(int n, float x0, float dx, int nwedge, float offset, bool inten)
3007{
3008 ENTERFUNC;
3009
3010 if (nz > 1) {
3011 LOGERR("2D images only.");
3012 throw ImageDimensionException("2D images only");
3013 }
3014 int isinten=get_attr_default("is_intensity",0);
3015
3016 if (isinten&&!inten) { throw InvalidParameterException("Must set inten for calc_radial_dist with intensity image"); }
3017
3018
3019 vector<double>ret(n*nwedge);
3020 vector<double>norm(n*nwedge);
3021
3022 int x,y,i;
3023 int isri = is_ri(); // this has become expensive!
3024 int step=is_complex()?2:1;
3025 float astep=(float)(M_PI*2.0/nwedge);
3026 if (is_complex()) astep/=2; // Since we only have the right 1/2 of Fourier space
3027 float* data = get_data();
3028 for (i=0; i<n*nwedge; i++) ret[i]=norm[i]=0.0;
3029
3030 // We do 2D separately to avoid the hypot3 call
3031 for (y=i=0; y<ny; y++) {
3032 for (x=0; x<nx; x+=step,i+=step) {
3033 float r,v,a;
3034 int bin;
3035 if (is_complex()) {
3036 r=(float)(hypot(x/2.0,y<ny/2?y:ny-y)); // origin at 0,0; periodic
3037 a=atan2(float(y<ny/2?y:y-ny),x/2.0f);
3038 if (!inten) {
3039 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
3040 else v=data[i]; // amp/phase, just get amp
3041 } else {
3042 if (isinten) v=data[i];
3043 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
3044 else v=data[i]*data[i];
3045 }
3046 bin=n*int(floor((a+M_PI/2.0f+offset)/astep));
3047 }
3048 else {
3049 r=(float)(hypot(x-nx/2,y-ny/2));
3050 a=atan2(float(y-ny/2),float(x-nx/2));
3051 if (inten) v=data[i]*data[i];
3052 else v=data[i];
3053 bin=n*int(floor((a+M_PI+offset)/astep));
3054 }
3055 if (bin>=nwedge*n) bin-=nwedge*n;
3056 if (bin<0) bin+=nwedge*n;
3057 r=(r-x0)/dx;
3058 int f=int(r); // safe truncation, so floor isn't needed
3059 r-=float(f); // r is now the fractional spacing between bins
3060// printf("%d %d %d %d %1.3f %1.3f\n",x,y,bin,f,r,a);
3061 if (f>=0 && f<n) {
3062 ret[f+bin]+=v*(1.0f-r);
3063 norm[f+bin]+=(1.0f-r);
3064 if (f<n-1) {
3065 ret[f+1+bin]+=v*r;
3066 norm[f+1+bin]+=r;
3067 }
3068 }
3069 }
3070 }
3071
3072 for (i=0; i<n*nwedge; i++) ret[i]/=norm[i]?norm[i]:1.0f; // Normalize
3073 EXITFUNC;
3074
3075 return vector<float>(ret.begin(),ret.end());
3076}
3077
3078void EMData::cconj() {
3079 ENTERFUNC;
3080 if (!is_complex() || !is_ri())
3081 throw ImageFormatException("EMData::conj requires a complex, ri image");
3082 int nxreal = nx -2 + int(is_fftodd());
3083 int nxhalf = nxreal/2;
3084 for (int iz = 0; iz < nz; iz++)
3085 for (int iy = 0; iy < ny; iy++)
3086 for (int ix = 0; ix <= nxhalf; ix++)
3087 cmplx(ix,iy,iz) = conj(cmplx(ix,iy,iz));
3088 EXITFUNC;
3089}
3090
3091void EMData::update_stat() const
3092{
3093 ENTERFUNC;
3094// printf("update stat %f %d\n",(float)attr_dict["mean"],flags);
3095 if (!(flags & EMDATA_NEEDUPD))
3096 {
3097 EXITFUNC;
3098 return;
3099 }
3100 if (rdata==0) return;
3101
3102 float* data = get_data();
3103 float max = -FLT_MAX;
3104 float min = -max;
3105
3106 double sum = 0;
3107 double square_sum = 0;
3108
3109 int step = 1;
3110 if (is_complex() && !is_ri()) {
3111 step = 2;
3112 }
3113
3114 int n_nonzero = 0;
3115
3116 size_t size = (size_t)nx*ny*nz;
3117
3118 for (size_t i = 0; i < size; i += step) {
3119 float v = data[i];
3120 max = Util::get_max(max, v);
3121 min = Util::get_min(min, v);
3122 sum += v;
3123 square_sum += v * (double)(v);
3124 if (v != 0) n_nonzero++;
3125 }
3126
3127 size_t n = size / step;
3128 double mean = sum / n;
3129 double var = (square_sum - sum*sum / n) / (n-1);
3130 double sigma = var >= 0.0 ? std::sqrt(var) : 0.0;
3131 n_nonzero = Util::get_max(1, n_nonzero);
3132 double varn = (square_sum - sum*sum / n_nonzero) / (n_nonzero-1);
3133 double sigma_nonzero = varn >= 0.0 ? std::sqrt(varn) : 0.0;
3134 double mean_nonzero = sum / n_nonzero; // previous version overcounted! G2
3135
3136 attr_dict["minimum"] = min;
3137 attr_dict["maximum"] = max;
3138 attr_dict["mean"] = (float)(mean);
3139 attr_dict["sigma"] = (float)(sigma);
3140 attr_dict["square_sum"] = (float)(square_sum);
3141 attr_dict["mean_nonzero"] = (float)(mean_nonzero);
3142 attr_dict["sigma_nonzero"] = (float)(sigma_nonzero);
3143 attr_dict["is_complex"] = (int) is_complex();
3144 attr_dict["is_complex_ri"] = (int) is_ri();
3145
3146 flags &= ~EMDATA_NEEDUPD;
3147
3148 if (rot_fp != 0)
3149 {
3150 delete rot_fp; rot_fp = 0;
3151 }
3152
3153 EXITFUNC;
3154// printf("done stat %f %f %f\n",(float)mean,(float)max,(float)sigma);
3155}
3156
3161bool EMData::operator==(const EMData& that) const {
3162 if(this != &that) {
3163 return false;
3164 }
3165 else {
3166 return true;
3167 }
3168}
3169
3170bool EMData::equal(const EMData& that) const {
3171 if (that.get_xsize() != nx || that.get_ysize() != ny || that.get_zsize() != nz ) return false;
3172
3173 const float* d1 = that.get_const_data();
3174 float* d2 = get_data();
3175
3176 for(size_t i =0; i < get_size(); ++i,++d1,++d2) {
3177 if ((*d1) != (*d2)) return false;
3178 }
3179
3180// if(attr_dict != that.attr_dict) {
3181// return false;
3182// }
3183
3184 return true;
3185}
3186
3187EMData * EMAN::operator+(const EMData & em, float n)
3188{
3189 EMData * r = em.copy();
3190 r->add(n);
3191 return r;
3192}
3193
3194EMData * EMAN::operator-(const EMData & em, float n)
3195{
3196 EMData* r = em.copy();
3197 r->sub(n);
3198 return r;
3199}
3200
3201EMData * EMAN::operator*(const EMData & em, float n)
3202{
3203 EMData* r = em.copy();
3204 r ->mult(n);
3205 return r;
3206}
3207
3208EMData * EMAN::operator/(const EMData & em, float n)
3209{
3210 EMData * r = em.copy();
3211 r->div(n);
3212 return r;
3213}
3214
3215
3216EMData * EMAN::operator+(float n, const EMData & em)
3217{
3218 EMData * r = em.copy();
3219 r->add(n);
3220 return r;
3221}
3222
3223EMData * EMAN::operator-(float n, const EMData & em)
3224{
3225 EMData * r = em.copy();
3226 r->mult(-1.0f);
3227 r->add(n);
3228 return r;
3229}
3230
3231EMData * EMAN::operator*(float n, const EMData & em)
3232{
3233 EMData * r = em.copy();
3234 r->mult(n);
3235 return r;
3236}
3237
3238EMData * EMAN::operator/(float n, const EMData & em)
3239{
3240 EMData * r = em.copy();
3241 r->to_one();
3242 r->mult(n);
3243 r->div(em);
3244
3245 return r;
3246}
3247
3248EMData * EMAN::rsub(const EMData & em, float n)
3249{
3250 return EMAN::operator-(n, em);
3251}
3252
3253EMData * EMAN::rdiv(const EMData & em, float n)
3254{
3255 return EMAN::operator/(n, em);
3256}
3257
3258EMData * EMAN::operator+(const EMData & a, const EMData & b)
3259{
3260 EMData * r = a.copy();
3261 r->add(b);
3262 return r;
3263}
3264
3265EMData * EMAN::operator-(const EMData & a, const EMData & b)
3266{
3267 EMData * r = a.copy();
3268 r->sub(b);
3269 return r;
3270}
3271
3272EMData * EMAN::operator*(const EMData & a, const EMData & b)
3273{
3274 EMData * r = a.copy();
3275 r->mult(b);
3276 return r;
3277}
3278
3279EMData * EMAN::operator/(const EMData & a, const EMData & b)
3280{
3281 EMData * r = a.copy();
3282 r->div(b);
3283 return r;
3284}
3285
3286void EMData::set_xyz_origin(float origin_x, float origin_y, float origin_z)
3287{
3288 attr_dict["origin_x"] = origin_x;
3289 attr_dict["origin_y"] = origin_y;
3290 attr_dict["origin_z"] = origin_z;
3291}
3292
3293#if 0
3294void EMData::calc_rcf(EMData * with, vector < float >&sum_array)
3295{
3296 ENTERFUNC;
3297
3298 int array_size = sum_array.size();
3299 float da = 2 * M_PI / array_size;
3300 float *dat = new float[array_size + 2];
3301 float *dat2 = new float[array_size + 2];
3302 int nx2 = nx * 9 / 20;
3303
3304 float lim = 0;
3305 if (fabs(translation[0]) < fabs(translation[1])) {
3306 lim = fabs(translation[1]);
3307 }
3308 else {
3309 lim = fabs(translation[0]);
3310 }
3311
3312 nx2 -= (int)(floor(lim));
3313
3314 for (int i = 0; i < array_size; i++) {
3315 sum_array[i] = 0;
3316 }
3317
3318 float sigma = attr_dict["sigma"];
3319 float with_sigma = with->get_attr_dict().get("sigma");
3320
3321 vector<float> vdata, vdata2;
3322 for (int i = 8; i < nx2; i += 6) {
3323 vdata = calc_az_dist(array_size, 0, da, i, i + 6);
3324 vdata2 = with->calc_az_dist(array_size, 0, da, i, i + 6);
3325 Assert(vdata.size() <= array_size + 2);
3326 Assert(cdata2.size() <= array_size + 2);
3327 std::copy(vdata.begin(), vdata.end(), dat);
3328 std::copy(vdata2.begin(), vdata2.end(), dat2);
3329
3330 EMfft::real_to_complex_1d(dat, dat, array_size);
3331 EMfft::real_to_complex_1d(dat2, dat2, array_size);
3332
3333 for (int j = 0; j < array_size + 2; j += 2) {
3334 float max = dat[j] * dat2[j] + dat[j + 1] * dat2[j + 1];
3335 float max2 = dat[j + 1] * dat2[j] - dat2[j + 1] * dat[j];
3336 dat[j] = max;
3337 dat[j + 1] = max2;
3338 }
3339
3340 EMfft::complex_to_real_1d(dat, dat, array_size);
3341 float norm = array_size * array_size * (4.0f * sigma) * (4.0f * with_sigma);
3342
3343 for (int j = 0; j < array_size; j++) {
3344 sum_array[j] += dat[j] * (float) i / norm;
3345 }
3346 }
3347
3348 if( dat )
3349 {
3350 delete[]dat;
3351 dat = 0;
3352 }
3353
3354 if( dat2 )
3355 {
3356 delete[]dat2;
3357 dat2 = 0;
3358 }
3359 EXITFUNC;
3360}
3361
3362#endif
3363
3364void EMData::add_incoherent(EMData * obj)
3365{
3366 ENTERFUNC;
3367
3368 if (!obj) {
3369 LOGERR("NULL image");
3370 throw NullPointerException("NULL image");
3371 }
3372
3373 if (!obj->is_complex() || !is_complex()) {
3374 throw ImageFormatException("complex images only");
3375 }
3376
3377 if (!EMUtil::is_same_size(this, obj)) {
3378 throw ImageFormatException("images not same size");
3379 }
3380
3381 ri2ap();
3382 obj->ri2ap();
3383
3384 float *dest = get_data();
3385 float *src = obj->get_data();
3386 size_t size = (size_t)nx * ny * nz;
3387 for (size_t j = 0; j < size; j += 2) {
3388 dest[j] = (float) hypot(src[j], dest[j]);
3389 dest[j + 1] = 0;
3390 }
3391
3392 obj->update();
3393 update();
3394 EXITFUNC;
3395}
3396
3397
3398float EMData::calc_dist(EMData * second_img, int y_index) const
3399{
3400 ENTERFUNC;
3401
3402 if (get_ndim() != 1) {
3403 throw ImageDimensionException("'this' image is 1D only");
3404 }
3405
3406 if (second_img->get_xsize() != nx || ny != 1) {
3407 throw ImageFormatException("image xsize not same");
3408 }
3409
3410 if (y_index > second_img->get_ysize() || y_index < 0) {
3411 return -1;
3412 }
3413
3414 float ret = 0;
3415 float *d1 = get_data();
3416 float *d2 = second_img->get_data() + second_img->get_xsize() * y_index;
3417
3418 for (int i = 0; i < nx; i++) {
3419 ret += Util::square(d1[i] - d2[i]);
3420 }
3421 EXITFUNC;
3422 return std::sqrt(ret);
3423}
3424
3425
3426EMData * EMData::calc_fast_sigma_image( EMData* mask)
3427{
3428 ENTERFUNC;
3429
3430 bool maskflag = false;
3431 if (mask == 0) {
3432 mask = new EMData(nx,ny,nz);
3433 mask->process_inplace("testimage.circlesphere");
3434 maskflag = true;
3435 }
3436
3437 if (get_ndim() != mask->get_ndim() ) throw ImageDimensionException("The dimensions do not match");
3438
3439 int mnx = mask->get_xsize(); int mny = mask->get_ysize(); int mnz = mask->get_zsize();
3440
3441 if ( mnx > nx || mny > ny || mnz > nz)
3442 throw ImageDimensionException("Can not calculate variance map using an image that is larger than this image");
3443
3444 size_t P = 0;
3445 for(size_t i = 0; i < mask->get_size(); ++i){
3446 if (mask->get_value_at(i) != 0){
3447 ++P;
3448 }
3449 }
3450 float normfac = 1.0f/(float)P;
3451
3452// bool undoclip = false;
3453
3454 int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
3455// if ( mnx < nx || mny < ny || mnz < nz) {
3456 Region r;
3457 if (ny == 1) r = Region((mnx-nxc)/2,nxc);
3458 else if (nz == 1) r = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
3459 else r = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
3460 mask->clip_inplace(r,0.0);
3461 //Region r((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
3462 //mask->clip_inplace(r);
3463 //undoclip = true;
3464 //}
3465
3466 // Here we generate the local average of the squares
3467 Region r2;
3468 if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
3469 else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
3470 else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
3471 EMData* squared = get_clip(r2,get_edge_mean());
3472
3473 EMData* tmp = squared->copy();
3474 Dict pow;
3475 pow["pow"] = 2.0f;
3476 squared->process_inplace("math.pow",pow);
3477 EMData* s = mask->convolute(squared);//ming, mask squared exchange
3478 squared->mult(normfac);
3479
3480 EMData* m = mask->convolute(tmp);//ming, tmp mask exchange
3481 m->mult(normfac);
3482 m->process_inplace("math.pow",pow);
3483 delete tmp; tmp = 0;
3484 s->sub(*m);
3485 // Here we finally generate the standard deviation image
3486 s->process_inplace("math.sqrt");
3487
3488// if ( undoclip ) {
3489// Region r((nx-mnx)/2, (ny-mny)/2, (nz-mnz)/2,mnx,mny,mnz);
3490// mask->clip_inplace(r);
3491// }
3492
3493 if (maskflag) {
3494 delete mask;
3495 mask = 0;
3496 } else {
3497 Region r;
3498 if (ny == 1) r = Region((nxc-mnx)/2,mnx);
3499 else if (nz == 1) r = Region((nxc-mnx)/2, (nyc-mny)/2,mnx,mny);
3500 else r = Region((nxc-mnx)/2, (nyc-mny)/2,(nzc-mnz)/2,mnx,mny,mnz);
3501 mask->clip_inplace(r);
3502 }
3503
3504 delete squared;
3505 delete m;
3506
3507 s->process_inplace("xform.phaseorigin.tocenter");
3508 Region r3;
3509 if (ny == 1) r3 = Region((nxc-nx)/2,nx);
3510 else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
3511 else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
3512 s->clip_inplace(r3);
3513 EXITFUNC;
3514 return s;
3515}
3516
3517// The following code looks strange - does anybody know it? Please let me know, pawel.a.penczek@uth.tmc.edu 04/09/06.
3518// This is just an implementation of "Roseman's" fast normalized cross-correlation (Ultramicroscopy, 2003). But the contents of this function have changed dramatically since you wrote that comment (d.woolford).
3519EMData *EMData::calc_flcf(EMData * with)
3520{
3521 ENTERFUNC;
3522 EMData *this_copy=this;
3523 this_copy=copy();
3524
3525 int mnx = with->get_xsize(); int mny = with->get_ysize(); int mnz = with->get_zsize();
3526 int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
3527
3528 // Ones is a circular/spherical mask, consisting of 1s.
3529 EMData* ones = new EMData(mnx,mny,mnz);
3530 ones->process_inplace("testimage.circlesphere");
3531
3532 // Get a copy of with, we will eventually resize it
3533 EMData* with_resized = with->copy();
3534 with_resized->process_inplace("normalize");
3535 with_resized->mult(*ones);
3536
3537 EMData* s = calc_fast_sigma_image(ones);// Get the local sigma image
3538
3539 Region r1;
3540 if (ny == 1) r1 = Region((mnx-nxc)/2,nxc);
3541 else if (nz == 1) r1 = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
3542 else r1 = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
3543 with_resized->clip_inplace(r1,0.0);
3544
3545 Region r2;
3546 if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
3547 else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
3548 else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
3549 this_copy->clip_inplace(r2,0.0);
3550
3551 EMData* corr = this_copy->calc_ccf(with_resized); // the ccf results should have same size as sigma
3552
3553 corr->process_inplace("xform.phaseorigin.tocenter");
3554 Region r3;
3555 if (ny == 1) r3 = Region((nxc-nx)/2,nx);
3556 else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
3557 else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
3558 corr->clip_inplace(r3);
3559
3560 corr->div(*s);
3561
3562 delete with_resized; delete ones; delete this_copy; delete s;
3563 EXITFUNC;
3564 return corr;
3565}
3566
3567EMData *EMData::convolute(EMData * with)
3568{
3569 ENTERFUNC;
3570
3571 EMData *f1 = do_fft();
3572 if (!f1) {
3573 LOGERR("FFT returns NULL image");
3574 throw NullPointerException("FFT returns NULL image");
3575 }
3576
3577 f1->ap2ri();
3578
3579 EMData *cf = 0;
3580 if (with) {
3581 if (with-> is_real())
3582 cf = with->do_fft();
3583 else
3584 cf = with->copy();
3585 if (!cf) {
3586 LOGERR("FFT returns NULL image");
3587 throw NullPointerException("FFT returns NULL image");
3588 }
3589 cf->ap2ri();
3590 }
3591 else {
3592 cf = f1->copy();
3593 }
3594 //printf("cf_x=%d, f1y=%d, thisx=%d, withx=%d\n",cf->get_xsize(),f1->get_ysize(),this->get_xsize(),with->get_xsize());
3595 if (with && !EMUtil::is_same_size(f1, cf)) {
3596 LOGERR("images not same size");
3597 throw ImageFormatException("images not same size");
3598 }
3599
3600 float *rdata1 = f1->get_data();
3601 float *rdata2 = cf->get_data();
3602 size_t cf_size = (size_t)cf->get_xsize() * cf->get_ysize() * cf->get_zsize();
3603
3604 float re,im;
3605
3606 for (size_t i = 0; i < cf_size; i += 2) {
3607 re = rdata1[i] * rdata2[i] - rdata1[i + 1] * rdata2[i + 1];
3608 im = rdata1[i + 1] * rdata2[i] + rdata1[i] * rdata2[i + 1];
3609 rdata2[i]=re;
3610 rdata2[i+1]=im;
3611 }
3612 cf->update();
3613 EMData *f2 = cf->do_ift();//ming change cf to cf_temp
3614 //printf("cf_x=%d, f2x=%d, thisx=%d, withx=%d\n",cf->get_xsize(),f2->get_xsize(),this->get_xsize(),with->get_xsize());
3615 if( cf )
3616 {
3617 delete cf;
3618 cf = 0;
3619 }
3620
3621 if( f1 )
3622 {
3623 delete f1;
3624 f1=0;
3625 }
3626
3627 EXITFUNC;
3628 return f2;
3629}
3630
3631
3632void EMData::common_lines(EMData * image1, EMData * image2,
3633 int mode, int steps, bool horizontal)
3634{
3635 ENTERFUNC;
3636
3637 if (!image1 || !image2) {
3638 throw NullPointerException("NULL image");
3639 }
3640
3641 if (mode < 0 || mode > 2) {
3642 throw OutofRangeException(0, 2, mode, "invalid mode");
3643 }
3644
3645 if (!image1->is_complex()) {
3646 image1 = image1->do_fft();
3647 }
3648 if (!image2->is_complex()) {
3649 image2 = image2->do_fft();
3650 }
3651
3652 image1->ap2ri();
3653 image2->ap2ri();
3654
3655 if (!EMUtil::is_same_size(image1, image2)) {
3656 throw ImageFormatException("images not same sizes");
3657 }
3658
3659 int image2_nx = image2->get_xsize();
3660 int image2_ny = image2->get_ysize();
3661
3662 int rmax = image2_ny / 4 - 1;
3663 int array_size = steps * rmax * 2;
3664 float *im1 = new float[array_size];
3665 float *im2 = new float[array_size];
3666 for (int i = 0; i < array_size; i++) {
3667 im1[i] = 0;
3668 im2[i] = 0;
3669 }
3670
3671 set_size(steps * 2, steps * 2, 1);
3672
3673 float *image1_data = image1->get_data();
3674 float *image2_data = image2->get_data();
3675
3676 float da = M_PI / steps;
3677 float a = -M_PI / 2.0f + da / 2.0f;
3678 int jmax = 0;
3679
3680 for (int i = 0; i < steps * 2; i += 2, a += da) {
3681 float s1 = 0;
3682 float s2 = 0;
3683 int i2 = i * rmax;
3684 int j = 0;
3685
3686 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) {
3687 float x = r * cos(a);
3688 float y = r * sin(a);
3689
3690 if (x < 0) {
3691 x = -x;
3692 y = -y;
3693 LOGERR("CCL ERROR %d, %f !\n", i, -x);
3694 }
3695
3696 int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx);
3697 int l = i2 + j;
3698 float x2 = x - floor(x);
3699 float y2 = y - floor(y);
3700
3701 im1[l] = Util::bilinear_interpolate(image1_data[k],
3702 image1_data[k + 2],
3703 image1_data[k + image2_nx],
3704 image1_data[k + 2 + image2_nx], x2, y2);
3705
3706 im2[l] = Util::bilinear_interpolate(image2_data[k],
3707 image2_data[k + 2],
3708 image2_data[k + image2_nx],
3709 image2_data[k + 2 + image2_nx], x2, y2);
3710
3711 k++;
3712
3713 im1[l + 1] = Util::bilinear_interpolate(image1_data[k],
3714 image1_data[k + 2],
3715 image1_data[k + image2_nx],
3716 image1_data[k + 2 + image2_nx], x2, y2);
3717
3718 im2[l + 1] = Util::bilinear_interpolate(image2_data[k],
3719 image2_data[k + 2],
3720 image2_data[k + image2_nx],
3721 image2_data[k + 2 + image2_nx], x2, y2);
3722
3723 s1 += Util::square_sum(im1[l], im1[l + 1]);
3724 s2 += Util::square_sum(im2[l], im2[l + 1]);
3725 }
3726
3727 jmax = j - 1;
3728 float sqrt_s1 = std::sqrt(s1);
3729 float sqrt_s2 = std::sqrt(s2);
3730
3731 int l = 0;
3732 for (float r = 1; r < rmax; r += 1.0f) {
3733 int i3 = i2 + l;
3734 im1[i3] /= sqrt_s1;
3735 im1[i3 + 1] /= sqrt_s1;
3736 im2[i3] /= sqrt_s2;
3737 im2[i3 + 1] /= sqrt_s2;
3738 l += 2;
3739 }
3740 }
3741 float * data = get_data();
3742
3743 switch (mode) {
3744 case 0:
3745 for (int m1 = 0; m1 < 2; m1++) {
3746 for (int m2 = 0; m2 < 2; m2++) {
3747
3748 if (m1 == 0 && m2 == 0) {
3749 for (int i = 0; i < steps; i++) {
3750 int i2 = i * rmax * 2;
3751 for (int j = 0; j < steps; j++) {
3752 int l = i + j * steps * 2;
3753 int j2 = j * rmax * 2;
3754 data[l] = 0;
3755 for (int k = 0; k < jmax; k++) {
3756 data[l] += im1[i2 + k] * im2[j2 + k];
3757 }
3758 }
3759 }
3760 }
3761 else {
3762 int steps2 = steps * m2 + steps * steps * 2 * m1;
3763
3764 for (int i = 0; i < steps; i++) {
3765 int i2 = i * rmax * 2;
3766 for (int j = 0; j < steps; j++) {
3767 int j2 = j * rmax * 2;
3768 int l = i + j * steps * 2 + steps2;
3769 data[l] = 0;
3770
3771 for (int k = 0; k < jmax; k += 2) {
3772 i2 += k;
3773 j2 += k;
3774 data[l] += im1[i2] * im2[j2];
3775 data[l] += -im1[i2 + 1] * im2[j2 + 1];
3776 }
3777 }
3778 }
3779 }
3780 }
3781 }
3782
3783 break;
3784 case 1:
3785 for (int m1 = 0; m1 < 2; m1++) {
3786 for (int m2 = 0; m2 < 2; m2++) {
3787 int steps2 = steps * m2 + steps * steps * 2 * m1;
3788 int p1_sign = 1;
3789 if (m1 != m2) {
3790 p1_sign = -1;
3791 }
3792
3793 for (int i = 0; i < steps; i++) {
3794 int i2 = i * rmax * 2;
3795
3796 for (int j = 0; j < steps; j++) {
3797 int j2 = j * rmax * 2;
3798
3799 int l = i + j * steps * 2 + steps2;
3800 data[l] = 0;
3801 float a = 0;
3802
3803 for (int k = 0; k < jmax; k += 2) {
3804 i2 += k;
3805 j2 += k;
3806
3807 float a1 = (float) hypot(im1[i2], im1[i2 + 1]);
3808 float p1 = atan2(im1[i2 + 1], im1[i2]);
3809 float p2 = atan2(im2[j2 + 1], im2[j2]);
3810
3811 data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1;
3812 a += a1;
3813 }
3814
3815 data[l] /= (float)(a * M_PI / 180.0f);
3816 }
3817 }
3818 }
3819 }
3820
3821 break;
3822 case 2:
3823 for (int m1 = 0; m1 < 2; m1++) {
3824 for (int m2 = 0; m2 < 2; m2++) {
3825 int steps2 = steps * m2 + steps * steps * 2 * m1;
3826
3827 for (int i = 0; i < steps; i++) {
3828 int i2 = i * rmax * 2;
3829
3830 for (int j = 0; j < steps; j++) {
3831 int j2 = j * rmax * 2;
3832 int l = i + j * steps * 2 + steps2;
3833 data[l] = 0;
3834
3835 for (int k = 0; k < jmax; k += 2) {
3836 i2 += k;
3837 j2 += k;
3838 data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1]));
3839 }
3840 }
3841 }
3842 }
3843 }
3844
3845 break;
3846 default:
3847 break;
3848 }
3849
3850 if (horizontal) {
3851 float *tmp_array = new float[ny];
3852 for (int i = 1; i < nx; i++) {
3853 for (int j = 0; j < ny; j++) {
3854 tmp_array[j] = get_value_at(i, j);
3855 }
3856 for (int j = 0; j < ny; j++) {
3857 set_value_at(i, j, 0, tmp_array[(j + i) % ny]);
3858 }
3859 }
3860 if( tmp_array )
3861 {
3862 delete[]tmp_array;
3863 tmp_array = 0;
3864 }
3865 }
3866
3867 if( im1 )
3868 {
3869 delete[]im1;
3870 im1 = 0;
3871 }
3872
3873 if( im2 )
3874 {
3875 delete[] im2;
3876 im2 = 0;
3877 }
3878
3879
3880 image1->update();
3881 image2->update();
3882 if( image1 )
3883 {
3884 delete image1;
3885 image1 = 0;
3886 }
3887 if( image2 )
3888 {
3889 delete image2;
3890 image2 = 0;
3891 }
3892 update();
3893 EXITFUNC;
3894}
3895
3896
3897
3898void EMData::common_lines_real(EMData * image1, EMData * image2,
3899 int steps, bool horiz)
3900{
3901 ENTERFUNC;
3902
3903 if (!image1 || !image2) {
3904 throw NullPointerException("NULL image");
3905 }
3906
3907 if (!EMUtil::is_same_size(image1, image2)) {
3908 throw ImageFormatException("images not same size");
3909 }
3910
3911 int steps2 = steps * 2;
3912 int image_ny = image1->get_ysize();
3913 EMData *image1_copy = image1->copy();
3914 EMData *image2_copy = image2->copy();
3915
3916 float *im1 = new float[steps2 * image_ny];
3917 float *im2 = new float[steps2 * image_ny];
3918
3919 EMData *images[] = { image1_copy, image2_copy };
3920 float *ims[] = { im1, im2 };
3921
3922 for (int m = 0; m < 2; m++) {
3923 float *im = ims[m];
3924 float a = M_PI / steps2;
3925 Transform t(Dict("type","2d","alpha",-a));
3926 for (int i = 0; i < steps2; i++) {
3927 images[i]->transform(t);
3928 float *data = images[i]->get_data();
3929
3930 for (int j = 0; j < image_ny; j++) {
3931 float sum = 0;
3932 for (int k = 0; k < image_ny; k++) {
3933 sum += data[j * image_ny + k];
3934 }
3935 im[i * image_ny + j] = sum;
3936 }
3937
3938 float sum1 = 0;
3939 float sum2 = 0;
3940 for (int j = 0; j < image_ny; j++) {
3941 int l = i * image_ny + j;
3942 sum1 += im[l];
3943 sum2 += im[l] * im[l];
3944 }
3945
3946 float mean = sum1 / image_ny;
3947 float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1);
3948
3949 for (int j = 0; j < image_ny; j++) {
3950 int l = i * image_ny + j;
3951 im[l] = (im[l] - mean) / sigma;
3952 }
3953
3954 images[i]->update();
3955 a += M_PI / steps;
3956 }
3957 }
3958
3959 set_size(steps2, steps2, 1);
3960 float *data1 = get_data();
3961
3962 if (horiz) {
3963 for (int i = 0; i < steps2; i++) {
3964 for (int j = 0, l = i; j < steps2; j++, l++) {
3965 if (l == steps2) {
3966 l = 0;
3967 }
3968
3969 float sum = 0;
3970 for (int k = 0; k < image_ny; k++) {
3971 sum += im1[i * image_ny + k] * im2[l * image_ny + k];
3972 }
3973 data1[i + j * steps2] = sum;
3974 }
3975 }
3976 }
3977 else {
3978 for (int i = 0; i < steps2; i++) {
3979 for (int j = 0; j < steps2; j++) {
3980 float sum = 0;
3981 for (int k = 0; k < image_ny; k++) {
3982 sum += im1[i * image_ny + k] * im2[j * image_ny + k];
3983 }
3984 data1[i + j * steps2] = sum;
3985 }
3986 }
3987 }
3988
3989 update();
3990
3991 if( image1_copy )
3992 {
3993 delete image1_copy;
3994 image1_copy = 0;
3995 }
3996
3997 if( image2_copy )
3998 {
3999 delete image2_copy;
4000 image2_copy = 0;
4001 }
4002
4003 if( im1 )
4004 {
4005 delete[]im1;
4006 im1 = 0;
4007 }
4008
4009 if( im2 )
4010 {
4011 delete[]im2;
4012 im2 = 0;
4013 }
4014 EXITFUNC;
4015}
4016
4017
4018void EMData::cut_slice(const EMData *const map, const Transform& transform, bool interpolate)
4019{
4020 ENTERFUNC;
4021
4022 if (!map) throw NullPointerException("NULL image");
4023 // These restrictions should be ultimately restricted so that all that matters is get_ndim() = (map->get_ndim() -1)
4024 if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
4025 if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
4026 // Now check for complex images - this is really just being thorough
4027 if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
4028 if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
4029
4030
4031 float *sdata = map->get_data();
4032 float *ddata = get_data();
4033
4034 int map_nx = map->get_xsize();
4035 int map_ny = map->get_ysize();
4036 int map_nz = map->get_zsize();
4037 int map_nxy = map_nx * map_ny;
4038
4039 int ymax = ny/2;
4040 if ( ny % 2 == 1 ) ymax += 1;
4041 int xmax = nx/2;
4042 if ( nx % 2 == 1 ) xmax += 1;
4043 for (int y = -ny/2; y < ymax; y++) {
4044 for (int x = -nx/2; x < xmax; x++) {
4045 Vec3f coord(x,y,0);
4046 Vec3f soln = transform*coord;
4047
4048// float xx = (x+pretrans[0]) * (*ort)[0][0] + (y+pretrans[1]) * (*ort)[0][1] + pretrans[2] * (*ort)[0][2] + posttrans[0];
4049// float yy = (x+pretrans[0]) * (*ort)[1][0] + (y+pretrans[1]) * (*ort)[1][1] + pretrans[2] * (*ort)[1][2] + posttrans[1];
4050// float zz = (x+pretrans[0]) * (*ort)[2][0] + (y+pretrans[1]) * (*ort)[2][1] + pretrans[2] * (*ort)[2][2] + posttrans[2];
4051
4052
4053// xx += map_nx/2;
4054// yy += map_ny/2;
4055// zz += map_nz/2;
4056
4057 float xx = soln[0]+map_nx/2;
4058 float yy = soln[1]+map_ny/2;
4059 float zz = soln[2]+map_nz/2;
4060
4061 int l = (x+nx/2) + (y+ny/2) * nx;
4062
4063 float t = xx - floor(xx);
4064 float u = yy - floor(yy);
4065 float v = zz - floor(zz);
4066
4067 if (xx < 0 || yy < 0 || zz < 0 ) {
4068 ddata[l] = 0;
4069 continue;
4070 }
4071 if (interpolate) {
4072 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
4073 ddata[l] = 0;
4074 continue;
4075 }
4076 int k = (int) (Util::fast_floor(xx) + Util::fast_floor(yy) * map_nx + Util::fast_floor(zz) * map_nxy);
4077
4078
4079 if (xx < (map_nx - 1) && yy < (map_ny - 1) && zz < (map_nz - 1)) {
4080 ddata[l] = Util::trilinear_interpolate(sdata[k],
4081 sdata[k + 1], sdata[k + map_nx],sdata[k + map_nx + 1],
4082 sdata[k + map_nxy], sdata[k + map_nxy + 1], sdata[k + map_nx + map_nxy],
4083 sdata[k + map_nx + map_nxy + 1],t, u, v);
4084 }
4085 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) && zz == (map_nz - 1) ) {
4086 ddata[l] += sdata[k];
4087 }
4088 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) ) {
4089 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + map_nxy],v);
4090 }
4091 else if ( xx == (map_nx - 1) && zz == (map_nz - 1) ) {
4092 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + map_nx],u);
4093 }
4094 else if ( yy == (map_ny - 1) && zz == (map_nz - 1) ) {
4095 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + 1],t);
4096 }
4097 else if ( xx == (map_nx - 1) ) {
4098 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + map_nx], sdata[k + map_nxy], sdata[k + map_nxy + map_nx],u,v);
4099 }
4100 else if ( yy == (map_ny - 1) ) {
4101 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nxy], sdata[k + map_nxy + 1],t,v);
4102 }
4103 else if ( zz == (map_nz - 1) ) {
4104 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nx], sdata[k + map_nx + 1],t,u);
4105 }
4106
4107// if (k >= map->get_size()) {
4108// cout << xx << " " << yy << " " << zz << " " << endl;
4109// cout << k << " " << get_size() << endl;
4110// cout << get_xsize() << " " << get_ysize() << " " << get_zsize() << endl;
4111// throw;
4112// }
4113//
4114// ddata[l] = Util::trilinear_interpolate(sdata[k],
4115// sdata[k + 1], sdata[k + map_nx],sdata[k + map_nx + 1],
4116// sdata[k + map_nxy], sdata[k + map_nxy + 1], sdata[k + map_nx + map_nxy],
4117// sdata[k + map_nx + map_nxy + 1],t, u, v);
4118 }
4119 else {
4120 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
4121 ddata[l] = 0;
4122 continue;
4123 }
4124 size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
4125 ddata[l] = sdata[k];
4126 }
4127
4128 }
4129 }
4130
4131 update();
4132
4133 EXITFUNC;
4134}
4135
4136EMData *EMData::unwrap_largerR(int r1,int r2,int xs, float rmax_f) {
4137 float *d,*dd;
4138 int do360=2;
4139 int rmax = (int)(rmax_f+0.5f);
4140 unsigned long i;
4141 unsigned int nvox=get_xsize()*get_ysize();//ming
4142 float maxmap=0.0f, minmap=0.0f;
4143 float temp=0.0f, diff_den=0.0f, norm=0.0f;
4144 float cut_off_va =0.0f;
4145
4146 d=get_data();
4147 maxmap=-1000000.0f;
4148 minmap=1000000.0f;
4149 for (i=0;i<nvox;i++){
4150 if(d[i]>maxmap) maxmap=d[i];
4151 if(d[i]<minmap) minmap=d[i];
4152 }
4153 diff_den = maxmap-minmap;
4154 for(i=0;i<nvox;i++) {
4155 temp = (d[i]-minmap)/diff_den;
4156 if(cut_off_va != 0.0) { // cut off the lowerset ?% noisy information
4157 if(temp < cut_off_va)
4158 d[i] = 0.0f; // set the empty part density=0.0
4159 else
4160 d[i] = temp-cut_off_va;
4161 }
4162 else d[i] = temp;
4163 }
4164
4165 for(i=0;i<nvox;i++) {
4166 temp=d[i];
4167 norm += temp*temp;
4168 }
4169 for(i=0;i<nvox;i++) d[i] /= norm; // y' = y/norm(y)
4170
4171 if (xs<1) {
4172 xs = (int) floor(do360*M_PI*get_ysize()/4); // ming
4173 xs=Util::calc_best_fft_size(xs); // ming
4174 }
4175 if (r1<0) r1=0;
4176 float maxext=ceil(0.6f*std::sqrt((float)(get_xsize()*get_xsize()+get_ysize()*get_ysize())));// ming add std::
4177
4178 if (r2<r1) r2=(int)maxext;
4179 EMData *ret = new EMData;
4180
4181 ret->set_size(xs,r2+1,1);
4182
4183 dd=ret->get_data();
4184
4185 for (int i=0; i<xs; i++) {
4186 float si=sin(i*M_PI*2/xs);
4187 float co=cos(i*M_PI*2/xs);
4188 for (int r=0; r<=maxext; r++) {
4189 float x=(r+r1)*co+get_xsize()/2; // ming
4190 float y=(r+r1)*si+get_ysize()/2; // ming
4191 if(x<0.0 || x>=get_xsize()-1.0 || y<0.0 || y>=get_ysize()-1.0 || r>rmax){ //Ming , ~~~~ rmax need pass here
4192 for(;r<=r2;r++) // here r2=MAXR
4193 dd[i+r*xs]=0.0;
4194 break;
4195 }
4196 int x1=(int)floor(x);
4197 int y1=(int)floor(y);
4198 float t=x-x1;
4199 float u=y-y1;
4200 float f11= d[x1+y1*get_xsize()]; // ming
4201 float f21= d[(x1+1)+y1*get_xsize()]; // ming
4202 float f12= d[x1+(y1+1)*get_xsize()]; // ming
4203 float f22= d[(x1+1)+(y1+1)*get_xsize()]; // ming
4204 dd[i+r*xs] = (1-t)*(1-u)*f11+t*(1-u)*f21+t*u*f22+(1-t)*u*f12;
4205 }
4206 }
4207 update();
4208 ret->update();
4209 return ret;
4210}
4211
4212
4213EMData *EMData::oneDfftPolar(int size, float rmax, float MAXR){ // sent MAXR value here later!!
4214 float *pcs=get_data();
4215 EMData *imagepcsfft = new EMData;
4216 imagepcsfft->set_size((size+2), (int)MAXR+1, 1);
4217 float *d=imagepcsfft->get_data();
4218
4219 EMData *data_in=new EMData;
4220 data_in->set_size(size,1,1);
4221 float *in=data_in->get_data();
4222
4223 for(int row=0; row<=(int)MAXR; ++row){
4224 if(row<=(int)rmax) {
4225 for(int i=0; i<size;++i) in[i] = pcs[i+row*size]; // ming
4226 data_in->set_complex(false);
4227 data_in->do_fft_inplace();
4228 for(int j=0;j<size+2;j++) d[j+row*(size+2)]=in[j];
4229 }
4230 else for(int j=0;j<size+2;j++) d[j+row*(size+2)]=0.0;
4231 }
4232 imagepcsfft->update();
4233 delete data_in;
4234 return imagepcsfft;
4235}
4236
4237void EMData::uncut_slice(EMData * const map, const Transform& transform) const
4238{
4239 ENTERFUNC;
4240
4241 if (!map) throw NullPointerException("NULL image");
4242 // These restrictions should be ultimately restricted so that all that matters is get_ndim() = (map->get_ndim() -1)
4243 if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
4244 if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
4245 // Now check for complex images - this is really just being thorough
4246 if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
4247 if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
4248
4249// Transform3D r( 0, 0, 0); // EMAN by default
4250// if (!ort) {
4251// ort = &r;
4252// }
4253
4254 float *ddata = map->get_data();
4255 float *sdata = get_data();
4256
4257 int map_nx = map->get_xsize();
4258 int map_ny = map->get_ysize();
4259 int map_nz = map->get_zsize();
4260 int map_nxy = map_nx * map_ny;
4261 float map_nz_round_limit = (float) map_nz-0.5f;
4262 float map_ny_round_limit = (float) map_ny-0.5f;
4263 float map_nx_round_limit = (float) map_nx-0.5f;
4264/*
4265 Vec3f posttrans = ort->get_posttrans();
4266 Vec3f pretrans = ort->get_pretrans();*/
4267
4268 int ymax = ny/2;
4269 if ( ny % 2 == 1 ) ymax += 1;
4270 int xmax = nx/2;
4271 if ( nx % 2 == 1 ) xmax += 1;
4272 for (int y = -ny/2; y < ymax; y++) {
4273 for (int x = -nx/2; x < xmax; x++) {
4274 Vec3f coord(x,y,0);
4275 Vec3f soln = transform*coord;
4276// float xx = (x+pretrans[0]) * (*ort)[0][0] + (y+pretrans[1]) * (*ort)[0][1] + pretrans[2] * (*ort)[0][2] + posttrans[0];
4277// float yy = (x+pretrans[0]) * (*ort)[1][0] + (y+pretrans[1]) * (*ort)[1][1] + pretrans[2] * (*ort)[1][2] + posttrans[1];
4278// float zz = (x+pretrans[0]) * (*ort)[2][0] + (y+pretrans[1]) * (*ort)[2][1] + pretrans[2] * (*ort)[2][2] + posttrans[2];
4279//
4280// xx += map_nx/2;
4281// yy += map_ny/2;
4282// zz += map_nz/2;
4283//
4284 float xx = soln[0]+map_nx/2;
4285 float yy = soln[1]+map_ny/2;
4286 float zz = soln[2]+map_nz/2;
4287
4288 // These 0.5 offsets are here because the round function rounds to the nearest whole number.
4289 if (xx < -0.5 || yy < -0.5 || zz < -0.5 || xx >= map_nx_round_limit || yy >= map_ny_round_limit || zz >= map_nz_round_limit) continue;
4290
4291 size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
4292 int l = (x+nx/2) + (y+ny/2) * nx;
4293 ddata[k] = sdata[l];
4294 }
4295 }
4296
4297 map->update();
4298 EXITFUNC;
4299}
4300
4301EMData *EMData::extract_box(const Transform& cs, const Region& r)
4302{
4303 vector<float> cs_matrix = cs.get_matrix();
4304
4305 EMData* box = new EMData();
4306 box->set_size((r.get_width()-r.x_origin()), (r.get_height()-r.y_origin()), (r.get_depth()-r.z_origin()));
4307 int box_nx = box->get_xsize();
4308 int box_ny = box->get_ysize();
4309 int box_nxy = box_nx*box_ny;
4310 float* bdata = box->get_data();
4311 float* ddata = get_data();
4312
4313 for (int x = r.x_origin(); x < r.get_width(); x++) {
4314 for (int y = r.y_origin(); y < r.get_height(); y++) {
4315 for (int z = r.z_origin(); z < r.get_depth(); z++) {
4316 //float xb = cs_matrix[0]*x + cs_matrix[1]*y + cs_matrix[2]*z + cs_matrix[3];
4317 //float yb = cs_matrix[4]*x + cs_matrix[5]*y + cs_matrix[6]*z + cs_matrix[7];
4318 //float zb = cs_matrix[8]*x + cs_matrix[9]*y + cs_matrix[10]*z + cs_matrix[11];
4319 float xb = cs_matrix[0]*x + y*cs_matrix[4] + z*cs_matrix[8] + cs_matrix[3];
4320 float yb = cs_matrix[1]*x + y*cs_matrix[5] + z*cs_matrix[9] + cs_matrix[7];
4321 float zb = cs_matrix[2]*x + y*cs_matrix[6] + z*cs_matrix[10] + cs_matrix[11];
4322 float t = xb - Util::fast_floor(xb);
4323 float u = yb - Util::fast_floor(yb);
4324 float v = zb - Util::fast_floor(zb);
4325
4326 //cout << x << " " << y << " " << z << " Box " << xb << " " << yb << " " << zb << endl;
4327 int l = (x - r.x_origin()) + (y - r.y_origin())*box_nx + (z - r.z_origin())*box_nxy;
4328 int k = (int) (Util::fast_floor(xb) + Util::fast_floor(yb) * nx + Util::fast_floor(zb) * nxy);
4329 //cout << k << " " << l << endl;
4330 if ( xb > nx - 1 || yb > ny - 1 || zb > nz - 1) {
4331 bdata[l] = 0;
4332 continue;
4333 }
4334 if (xb < 0 || yb < 0 || zb < 0){
4335 bdata[l] = 0;
4336 continue;
4337 }
4338
4339 if (xb < (nx - 1) && yb < (ny - 1) && zb < (nz - 1)) {
4340 bdata[l] = Util::trilinear_interpolate(ddata[k], ddata[k + 1], ddata[k + nx],ddata[k + nx + 1], ddata[k + nxy], ddata[k + nxy + 1], ddata[k + nx + nxy], ddata[k + nx + nxy + 1],t, u, v);
4341 }
4342 }
4343 }
4344 }
4345
4346 return box;
4347}
4348
4349void EMData::save_byteorder_to_dict(ImageIO * imageio)
4350{
4351 string image_endian = "ImageEndian";
4352 string host_endian = "HostEndian";
4353
4354 if (imageio->is_image_big_endian()) {
4355 attr_dict[image_endian] = "big";
4356 }
4357 else {
4358 attr_dict[image_endian] = "little";
4359 }
4360
4361 if (ByteOrder::is_host_big_endian()) {
4362 attr_dict[host_endian] = "big";
4363 }
4364 else {
4365 attr_dict[host_endian] = "little";
4366 }
4367}
4368
4369EMData* EMData::compute_missingwedge(float wedgeangle, float start, float stop)
4370{
4371 EMData* test = new EMData();
4372 test->set_size(nx,ny,nz);
4373
4374 float ratio = tan((90.0f-wedgeangle)*M_PI/180.0f);
4375
4376 int offset_i = 2*int(start*nz/2);
4377 int offset_f = int(stop*nz/2);
4378
4379 int step = 0;
4380 float sum = 0.0;
4381 double square_sum = 0.0;
4382 for (int j = 0; j < offset_f; j++){
4383 for (int k = offset_i; k < offset_f; k++) {
4384 for (int i = 0; i < nx; i+=2) {
4385 if (i < int(k*ratio)) {
4386 test->set_value_at(i, j, k, 1.0);
4387 float v = std::sqrt(pow(get_value_at_wrap(i, j, k),2) + pow(get_value_at_wrap(i+1, j, k),2));
4388 sum += v;
4389 square_sum += v * (double)(v);
4390 step++;
4391 }
4392 }
4393 }
4394 }
4395
4396 float mean = sum / step;
4397 double var = (square_sum - sum*mean) / (step-1);
4398 float sigma = var >= 0.0 ? (float) std::sqrt(var) : 0.0;
4399
4400 cout << "Mean sqr wedge amp " << mean << " Sigma Squ wedge Amp " << sigma << endl;
4401 set_attr("spt_wedge_mean", mean);
4402 set_attr("spt_wedge_sigma", sigma);
4403
4404 return test;
4405}
all_imageio.h
rdata
#define rdata(i)
Definition: analyzer.cpp:592
EMAN::ByteOrder::is_host_big_endian
static bool is_host_big_endian()
Definition: byteorder.cpp:40
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::Dict::has_key
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::EMData::get_top_half
EMData * get_top_half() const
Get the top half of this 3D image.
Definition: emdata.cpp:707
EMAN::EMData::scale
void scale(float scale_factor)
scale the image by a factor.
Definition: emdata.cpp:886
EMAN::EMData::update_stat
void update_stat() const
Definition: emdata.cpp:3091
EMAN::EMData::translate
void translate(float dx, float dy, float dz)
Translate this image.
Definition: emdata.cpp:904
EMAN::EMData::calc_fast_sigma_image
EMData * calc_fast_sigma_image(EMData *mask)
Calculates the local standard deviation (sigma) image using the given mask image.
Definition: emdata.cpp:3426
EMAN::EMData::pathnum
int pathnum
Definition: emdata.h:858
EMAN::EMData::compute_missingwedge
EMData * compute_missingwedge(float wedgeangle, float start=0.05, float stop=0.5)
Find the mean and variance of voxels in the missing wedge.
Definition: emdata.cpp:4369
EMAN::EMData::rot_fp
EMData * rot_fp
This is a cached rotational footprint, can save much time.
Definition: emdata.h:861
EMAN::EMData::cconj
void cconj()
Replace the image its complex conjugate.
Definition: emdata.cpp:3078
EMAN::EMData::rotate
void rotate(const Transform &t)
Rotate this image.
Definition: emdata.cpp:971
EMAN::EMData::calc_flcf
EMData * calc_flcf(EMData *with)
Calculates the cross correlation with local normalization between 2 images.
Definition: emdata.cpp:3519
EMAN::EMData::save_byteorder_to_dict
void save_byteorder_to_dict(ImageIO *imageio)
Definition: emdata.cpp:4349
EMAN::EMData::max_3D_pixel_error
float max_3D_pixel_error(const Transform &t1, const Transform &t2, float r)
Definition: emdata.cpp:977
EMAN::EMData::yoff
int yoff
Definition: emdata.h:851
EMAN::EMData::EMData
EMData()
For all image I/O.
Definition: emdata.cpp:70
EMAN::EMData::uncut_slice
void uncut_slice(EMData *const map, const Transform &tr) const
Opposite of the cut_slice().
Definition: emdata.cpp:4237
EMAN::EMData::window_center
EMData * window_center(int l)
Window the center of an image.
Definition: emdata.cpp:776
EMAN::EMData::totalalloc
static int totalalloc
Definition: emdata.h:803
EMAN::EMData::supp
float * supp
supplementary data array
Definition: emdata.h:837
EMAN::EMData::make_rotational_footprint_cmc
EMData * make_rotational_footprint_cmc(bool unwrap=true)
Definition: emdata.cpp:1795
EMAN::EMData::EMDATA_NEEDUPD
@ EMDATA_NEEDUPD
Definition: emdata.h:816
EMAN::EMData::calc_mutual_correlation
EMData * calc_mutual_correlation(EMData *with, bool tocorner=false, EMData *filter=0)
Calculates mutual correlation function (MCF) between 2 images.
Definition: emdata.cpp:2209
EMAN::EMData::changecount
int changecount
Definition: emdata.h:846
EMAN::EMData::zero_corner_circulant
void zero_corner_circulant(const int radius=0)
Zero the pixels in the bottom left corner of the image If radius is greater than 1,...
Definition: emdata.cpp:1462
EMAN::EMData::dot_rotate_translate
double dot_rotate_translate(EMData *with, float dx, float dy, float da, const bool mirror=false)
dot product of 2 images.
Definition: emdata.cpp:1239
EMAN::EMData::little_big_dot
EMData * little_big_dot(EMData *little_img, bool do_sigma=false)
This does a normalized dot product of a little image with a big image using real-space methods.
Definition: emdata.cpp:1322
EMAN::EMData::rdata
float * rdata
image real data
Definition: emdata.h:835
EMAN::EMData::calc_ccf_masked
EMData * calc_ccf_masked(EMData *with, EMData *withsq=0, EMData *mask=0)
Calculate cross correlation between this and with where this is assumed to have had a mask applied to...
Definition: emdata.cpp:1600
EMAN::EMData::attr_dict
Dict attr_dict
to store all image header info
Definition: emdata.h:833
EMAN::EMData::path
string path
Definition: emdata.h:857
EMAN::EMData::nx
int nx
image size
Definition: emdata.h:848
EMAN::EMData::rotate_x
void rotate_x(int dx)
This performs a translation of each line along x with wraparound.
Definition: emdata.cpp:1209
EMAN::EMData::get_clip
EMData * get_clip(const Region &area, const float fill=0) const
Get an inclusive clip.
Definition: emdata.cpp:592
EMAN::EMData::make_rotational_footprint_e1
EMData * make_rotational_footprint_e1(bool unwrap=true)
Definition: emdata.cpp:1898
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::apply_radial_func
void apply_radial_func(float x0, float dx, vector< float >array, bool interp=true)
multiplies by a radial function in fourier space.
Definition: emdata.cpp:2677
EMAN::EMData::calc_radial_dist
vector< float > calc_radial_dist(int n, float x0, float dx, int inten)
calculates radial distribution.
Definition: emdata.cpp:2781
EMAN::EMData::calc_hist
vector< float > calc_hist(int hist_size=128, float hist_min=0, float hist_max=0, const float &brt=0.0f, const float &cont=1.0f)
Calculates the histogram of 'this' image.
Definition: emdata.cpp:2365
EMAN::EMData::transform
void transform(const Transform &t)
Rotate then translate the image.
Definition: emdata.h:295
EMAN::EMData::cut_slice
void cut_slice(const EMData *const map, const Transform &tr, bool interpolate=true)
cut a 2D slice out of a real 3D map.
Definition: emdata.cpp:4018
EMAN::EMData::add_incoherent
void add_incoherent(EMData *obj)
Adds 'obj' to 'this' incoherently.
Definition: emdata.cpp:3364
EMAN::EMData::get_rotated_clip
EMData * get_rotated_clip(const Transform &xform, const IntSize &size, float scale=1.0)
This will extract an arbitrarily oriented and sized region from the image.
Definition: emdata.cpp:733
EMAN::EMData::do_radon
EMData * do_radon()
Radon Transform: an algorithm that transforms an original image into a series of equiangular projecti...
Definition: emdata.cpp:1413
EMAN::EMData::make_footprint
EMData * make_footprint(int type=0)
Makes a 'footprint' for the current image.
Definition: emdata.cpp:2035
EMAN::EMData::oneDfftPolar
EMData * oneDfftPolar(int size, float rmax, float MAXR)
Definition: emdata.cpp:4213
EMAN::EMData::all_translation
Vec3f all_translation
translation from the original location
Definition: emdata.h:854
EMAN::EMData::set_xyz_origin
void set_xyz_origin(float origin_x, float origin_y, float origin_z)
Set the x,y,z origin of the image.
Definition: emdata.cpp:3286
EMAN::EMData::convolute
EMData * convolute(EMData *with)
Convolutes 2 data sets.
Definition: emdata.cpp:3567
EMAN::EMData::calc_az_dist
vector< float > calc_az_dist(int n, float a0, float da, float rmin, float rmax)
Caculates the azimuthal distributions.
Definition: emdata.cpp:2472
EMAN::EMData::setup4slice
float * setup4slice(bool redo=true)
Set up for fftslice operations.
Definition: emdata.cpp:822
EMAN::EMData::unwrap_largerR
EMData * unwrap_largerR(int r1, int r2, int xs, float rmax_f)
Definition: emdata.cpp:4136
EMAN::EMData::extract_box
EMData * extract_box(const Transform &cs, const Region &r)
Extract a box from EMData in an abritrary orrientation.
Definition: emdata.cpp:4301
EMAN::EMData::rotate_translate
void rotate_translate(const Transform &t)
Apply a transformation to the image.
Definition: emdata.h:306
EMAN::EMData::clip_inplace
void clip_inplace(const Region &area, const float &fill_value=0)
Clip the image inplace - clipping region must be smaller than the current region internally memory is...
Definition: emdata.cpp:350
EMAN::EMData::calc_dist
float calc_dist(EMData *second_img, int y_index=0) const
Calculates the distance between 2 vectors.
Definition: emdata.cpp:3398
EMAN::EMData::xoff
int xoff
array index offsets
Definition: emdata.h:851
EMAN::EMData::zoff
int zoff
Definition: emdata.h:851
EMAN::EMUtil::em_memset
static void em_memset(void *data, const int value, const size_t size)
Definition: emutil.h:377
EMAN::EMUtil::em_memcpy
static void em_memcpy(void *dst, const void *const src, const size_t size)
Definition: emutil.h:384
EMAN::EMUtil::is_same_size
static bool is_same_size(const EMData *image1, const EMData *image2)
Check whether two EMData images are of the same size.
Definition: emutil.cpp:1224
EMAN::EMUtil::em_calloc
static void * em_calloc(const size_t nmemb, const size_t size)
Definition: emutil.h:370
EMAN::EMUtil::em_free
static void em_free(void *data)
Definition: emutil.h:380
EMAN::EMUtil::em_malloc
static void * em_malloc(const size_t size)
Definition: emutil.h:366
EMAN::ImageIO
ImageIO classes are designed for reading/writing various electron micrography image formats,...
Definition: imageio.h:127
EMAN::ImageIO::is_image_big_endian
virtual bool is_image_big_endian()=0
Is this image in big endian or not.
EMAN::IntPoint
IntPoint defines an integer-coordinate point in a 1D/2D/3D space.
Definition: geometry.h:192
EMAN::IntSize
IntSize is used to describe a 1D, 2D or 3D rectangular size in integers.
Definition: geometry.h:49
EMAN::Region
Region defines a 2D or 3D rectangular region specified by its origin coordinates and all edges' sizes...
Definition: geometry.h:497
EMAN::Region::get_width
float get_width() const
get the width
Definition: geometry.h:606
EMAN::Region::y_origin
float y_origin() const
get the y element of the origin
Definition: geometry.h:622
EMAN::Region::x_origin
float x_origin() const
get the x element of the origin
Definition: geometry.h:620
EMAN::Region::get_depth
float get_depth() const
get the depth
Definition: geometry.h:610
EMAN::Region::size
FloatSize size
Definition: geometry.h:655
EMAN::Region::origin
FloatPoint origin
Definition: geometry.h:654
EMAN::Region::z_origin
float z_origin() const
get the z element of the origin
Definition: geometry.h:624
EMAN::Region::get_height
float get_height() const
get the height
Definition: geometry.h:608
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::Transform::transform
Vec2f transform(const float &x, const float &y) const
Transform 2D coordinates using the internal transformation matrix.
Definition: transform.h:417
EMAN::Transform::set_rotation
void set_rotation(const Dict &rotation)
Set a rotation using a specific Euler type and the dictionary interface Works for all Euler types.
Definition: transform.cpp:519
EMAN::Transform::inverse
Transform inverse() const
Get the inverse of this transformation matrix.
Definition: transform.cpp:1327
EMAN::Transform::set_scale
void set_scale(const float &scale)
Set the scale.
Definition: transform.cpp:1123
EMAN::Transform::set_pre_trans
void set_pre_trans(const type &v)
Set the translational component of the matrix as though it was MSRT_ not MTSR, where T_ is the pre tr...
Definition: transform.h:556
EMAN::Transform::set_trans
void set_trans(const float &x, const float &y, const float &z=0)
Set the post translation component.
Definition: transform.cpp:1036
EMAN::Transform::get_matrix
vector< float > get_matrix() const
Get the transformation matrix using a vector.
Definition: transform.cpp:181
EMAN::Util::linear_interpolate
static float linear_interpolate(float p1, float p2, float t)
Calculate linear interpolation.
Definition: util.h:518
EMAN::Util::angle_sub_2pi
static float angle_sub_2pi(float x, float y)
Calculate the difference of 2 angles and makes the equivalent result to be less than Pi.
Definition: util.h:1066
EMAN::Util::square
static int square(int n)
Calculate a number's square.
Definition: util.h:736
EMAN::Util::fast_floor
static int fast_floor(float x)
A fast way to calculate a floor, which is largest integral value not greater than argument.
Definition: util.h:874
EMAN::Util::hypot3
static float hypot3(int x, int y, int z)
Euclidean distance function in 3D: f(x,y,z) = sqrt(x*x + y*y + z*z);.
Definition: util.h:827
EMAN::Util::round
static int round(float x)
Get ceiling round of a float number x.
Definition: util.h:465
EMAN::Util::calc_best_fft_size
static int calc_best_fft_size(int low)
Search the best FFT size with good primes.
Definition: util.cpp:1021
EMAN::Util::trilinear_interpolate
static float trilinear_interpolate(float p1, float p2, float p3, float p4, float p5, float p6, float p7, float p8, float t, float u, float v)
Calculate trilinear interpolation.
Definition: util.h:619
EMAN::Util::get_min
static int get_min(int f1, int f2)
Get the minimum of 2 numbers.
Definition: util.h:922
EMAN::Util::bilinear_interpolate
static float bilinear_interpolate(float p1, float p2, float p3, float p4, float t, float u)
Calculate bilinear interpolation.
Definition: util.h:543
EMAN::Util::square_sum
static float square_sum(float x, float y)
Calcuate (x*x + y*y).
Definition: util.h:764
EMAN::Util::get_max
static float get_max(float f1, float f2)
Get the maximum of 2 numbers.
Definition: util.h:998
EMAN::Util::hypot_fast
static float hypot_fast(int x, int y)
Euclidean distance in 2D for integers computed fast using a cached lookup table.
Definition: util.cpp:742
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)
cuda_processor.h
calc_ccf_cuda
void calc_ccf_cuda(float *afft, const float *bfft, const int nx, const int ny, const int nz)
mcf_cuda
void mcf_cuda(const float *data1, float *data2, const int nx, const int ny, const int nz)
emdata_unwrap
void emdata_unwrap(float *data, int r1, int r2, int xs, int num_pi, int dx, int dy, int weight_radial, int nx, int ny)
emdata_column_sum
float * emdata_column_sum(const float *data, const int nx, const int ny)
mult_complex_efficient_cuda
void mult_complex_efficient_cuda(float *data, const float *src_data, const int nx, const int ny, const int nz, const int radius)
calc_conv_cuda
void calc_conv_cuda(float *afft, const float *bfft, const int nx, const int ny, const int nz)
Assert
#define Assert(s)
Define Assert() function that is effective only when -DDEBUG is used.
Definition: emassert.h:42
sget_value_at_interp
float sget_value_at_interp(float x, float y) const
Get pixel density value at interpolation of (x,y).
Definition: emdata_core.cpp:914
copy
EMData * copy() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
to_zero
void to_zero()
Set all the pixel value = 0.
Definition: emdata_core.cpp:1465
cmplx
std::complex< float > & cmplx(const int ix, const int iy, const int iz)
Return reference to complex elements.
Definition: emdata_core.h:781
set_value_at
void set_value_at(Vec3i loc, float val)
set_value_at with Vec3i
Definition: emdata_core.h:552
equal
bool equal(const EMData &that) const
compare the equality of two EMData object based on their pixel values
get_value_at_wrap
float get_value_at_wrap(int x, int y, int z) const
Get the pixel density value at coordinates (x,y,z).
Definition: emdata_core.cpp:839
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
free_memory
void free_memory()
Free memory associated with this EMData Called in destructor and in assignment operator.
copy_head
EMData * copy_head() const
Make an image with a copy of the current image's header.
Definition: emdata_core.cpp:100
operator==
bool operator==(const EMData &that) const
read_image
void read_image(const string &filename, int img_index=0, bool header_only=false, const Region *region=0, bool is_3d=false, EMUtil::ImageType imgtype=EMUtil::IMAGE_UNKNOWN)
read an image file and stores its information to this EMData object.
is_fftodd
bool is_fftodd() const
Does this image correspond to a (real-space) odd nx?
Definition: emdata_metadata.h:864
get_attr
EMObject get_attr(const string &attr_name) const
The generic way to get any image header information given a header attribute name.
Definition: emdata_metadata.cpp:1006
get_const_data
const float * get_const_data() const
Get the image pixel density data in a 1D float array - const version of get_data.
Definition: emdata_metadata.h:84
is_fftpadded
bool is_fftpadded() const
Is this image already extended along x for ffts?
Definition: emdata_metadata.h:831
is_real
bool is_real() const
Is this a real image?
Definition: emdata_metadata.h:667
get_ysize
int get_ysize() const
Get the image y-dimensional size.
Definition: emdata_metadata.h:553
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_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
is_complex_x
bool is_complex_x() const
Is this image a 1D FFT image in X direction?
Definition: emdata_metadata.h:727
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_attr_default
EMObject get_attr_default(const string &attr_name, const EMObject &em_obj=EMObject()) const
The generic way to get any image header information given a header attribute name.
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.
get_size
size_t get_size() const
Get the number of allocated floats in the image (nx*ny*nz)
Definition: emdata_metadata.h:577
get_edge_mean
float get_edge_mean() const
Calculates the mean pixel values around the (1 pixel) edge of the image.
Definition: emdata_metadata.cpp:713
process
EMData * process(const string &processorname, const Dict &params=Dict()) const
Apply a processor with its parameters on a copy of this image, return result as a a new image.
process_inplace
void process_inplace(const string &processorname, const Dict &params=Dict())
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
ap2ri
void ap2ri()
convert the complex image from amplitude/phase to real/imaginary
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
InvalidParameterException
#define InvalidParameterException(desc)
Definition: exception.h:361
ImageFormatException
#define ImageFormatException(desc)
Definition: exception.h:147
ImageDimensionException
#define ImageDimensionException(desc)
Definition: exception.h:166
OutofRangeException
#define OutofRangeException(low, high, input, objname)
Definition: exception.h:334
UnexpectedBehaviorException
#define UnexpectedBehaviorException(desc)
Definition: exception.h:400
NullPointerException
#define NullPointerException(desc)
Definition: exception.h:241
EMAN::EMData::operator=
EMData & operator=(const EMData &that)
EMData assignment operator Performs a deep copy.
Definition: emdata.cpp:173
EMAN::EMData::calc_ccf
EMData * calc_ccf(EMData *with=0, fp_flag fpflag=CIRCULANT, bool center=false)
Calculate Cross-Correlation Function (CCF).
Definition: emdata.cpp:1499
EMAN::EMData::unwrap
EMData * unwrap(int r1=-1, int r2=-1, int xs=-1, int dx=0, int dy=0, bool do360=false, bool weight_radial=true) const
Maps to polar coordinates from Cartesian coordinates.
Definition: emdata.cpp:2596
EMAN::EMData::calc_ccfx
EMData * calc_ccfx(EMData *const with, int y0=0, int y1=-1, bool nosum=false, bool flip=false, bool usez=false)
Calculate Cross-Correlation Function (CCF) in the x-direction and adds them up, result in 1D.
Definition: emdata.cpp:1628
EMAN::EMData::make_rotational_footprint
EMData * make_rotational_footprint(bool unwrap=true)
Makes a 'rotational footprint', which is an 'unwound' autocorrelation function.
Definition: emdata.cpp:1862
EMAN::EMData::common_lines
void common_lines(EMData *image1, EMData *image2, int mode=0, int steps=180, bool horizontal=false)
Finds common lines between 2 complex images.
Definition: emdata.cpp:3632
EMAN::EMData::common_lines_real
void common_lines_real(EMData *image1, EMData *image2, int steps=180, bool horizontal=false)
Finds common lines between 2 real images.
Definition: emdata.cpp:3898
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
EMAN::operator+
EMData * operator+(const EMData &em, float n)
Definition: emdata.cpp:3187
EMAN::Vec3f
Vec3< float > Vec3f
Definition: vec3.h:693
EMAN::rdiv
EMData * rdiv(const EMData &em, float n)
Definition: emdata.cpp:3253
EMAN::Vec3i
Vec3< int > Vec3i
Definition: vec3.h:694
EMAN::dot
double dot(const Vector3 &w, const Vector3 &v)
Definition: vecmath.h:305
EMAN::rsub
EMData * rsub(const EMData &em, float n)
Definition: emdata.cpp:3248
EMAN::operator-
EMData * operator-(const EMData &em, float n)
Definition: emdata.cpp:3194
EMAN::operator/
EMData * operator/(const EMData &em, float n)
Definition: emdata.cpp:3208
EMAN::operator*
EMData * operator*(const EMData &em, float n)
Definition: emdata.cpp:3201
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517
images
#define images(i, j, k)
Definition: projector.cpp:1897
Generated by doxygen 1.9.2