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

code complete but contains bugs


Functions

void EMAN::EMData::common_lines (EMData *image1, EMData *image2, int mode=0, int steps=180, bool horizontal=false)
 Finds common lines between 2 complex images.
void EMAN::EMData::common_lines_real (EMData *image1, EMData *image2, int steps=180, bool horizontal=false)
 Finds common lines between 2 real images.

Function Documentation

void EMData::common_lines EMData image1,
EMData image2,
int  mode = 0,
int  steps = 180,
bool  horizontal = false
[inherited]
 

Finds common lines between 2 complex images.

This function does not assume any symmetry, just blindly compute the "sinogram" and the user has to take care how to interpret the returned "sinogram". it only considers inplane rotation and assumes prefect centering and identical scale.

Parameters:
image1 The first complex image.
image2 The second complex image.
mode Either 0 or 1 or 2. mode 0 is a summed dot-product, larger value means better match; mode 1 is weighted phase residual, lower value means better match.
steps 1/2 of the resolution of the map.
horizontal In horizontal way or not.
Exceptions:
NullPointerException If 'image1' or 'image2' is NULL.
OutofRangeException If 'mode' is invalid.
ImageFormatException If 'image1' 'image2' are not same size.

Definition at line 3442 of file emdata.cpp.

References EMAN::EMData::ap2ri(), data, EMAN::EMData::do_fft(), EMAN::EMData::get_data(), EMAN::EMData::get_value_at(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), ImageFormatException, EMAN::EMData::is_complex(), LOGERR, NullPointerException, ny, OutofRangeException, EMAN::EMData::set_size(), EMAN::EMData::set_value_at(), sqrt(), EMAN::EMData::update(), x, and y.

Referenced by main().

03444 {
03445         ENTERFUNC;
03446 
03447         if (!image1 || !image2) {
03448                 throw NullPointerException("NULL image");
03449         }
03450 
03451         if (mode < 0 || mode > 2) {
03452                 throw OutofRangeException(0, 2, mode, "invalid mode");
03453         }
03454 
03455         if (!image1->is_complex()) {
03456                 image1 = image1->do_fft();
03457         }
03458         if (!image2->is_complex()) {
03459                 image2 = image2->do_fft();
03460         }
03461 
03462         image1->ap2ri();
03463         image2->ap2ri();
03464 
03465         if (!EMUtil::is_same_size(image1, image2)) {
03466                 throw ImageFormatException("images not same sizes");
03467         }
03468 
03469         int image2_nx = image2->get_xsize();
03470         int image2_ny = image2->get_ysize();
03471 
03472         int rmax = image2_ny / 4 - 1;
03473         int array_size = steps * rmax * 2;
03474         float *im1 = new float[array_size];
03475         float *im2 = new float[array_size];
03476         for (int i = 0; i < array_size; i++) {
03477                 im1[i] = 0;
03478                 im2[i] = 0;
03479         }
03480 
03481         set_size(steps * 2, steps * 2, 1);
03482 
03483         float *image1_data = image1->get_data();
03484         float *image2_data = image2->get_data();
03485 
03486         float da = M_PI / steps;
03487         float a = -M_PI / 2.0f + da / 2.0f;
03488         int jmax = 0;
03489 
03490         for (int i = 0; i < steps * 2; i += 2, a += da) {
03491                 float s1 = 0;
03492                 float s2 = 0;
03493                 int i2 = i * rmax;
03494                 int j = 0;
03495 
03496                 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) {
03497                         float x = r * cos(a);
03498                         float y = r * sin(a);
03499 
03500                         if (x < 0) {
03501                                 x = -x;
03502                                 y = -y;
03503                                 LOGERR("CCL ERROR %d, %f !\n", i, -x);
03504                         }
03505 
03506                         int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx);
03507                         int l = i2 + j;
03508                         float x2 = x - floor(x);
03509                         float y2 = y - floor(y);
03510 
03511                         im1[l] = Util::bilinear_interpolate(image1_data[k],
03512                                                                                                 image1_data[k + 2],
03513                                                                                                 image1_data[k + image2_nx],
03514                                                                                                 image1_data[k + 2 + image2_nx], x2, y2);
03515 
03516                         im2[l] = Util::bilinear_interpolate(image2_data[k],
03517                                                                                                 image2_data[k + 2],
03518                                                                                                 image2_data[k + image2_nx],
03519                                                                                                 image2_data[k + 2 + image2_nx], x2, y2);
03520 
03521                         k++;
03522 
03523                         im1[l + 1] = Util::bilinear_interpolate(image1_data[k],
03524                                                                                                         image1_data[k + 2],
03525                                                                                                         image1_data[k + image2_nx],
03526                                                                                                         image1_data[k + 2 + image2_nx], x2, y2);
03527 
03528                         im2[l + 1] = Util::bilinear_interpolate(image2_data[k],
03529                                                                                                         image2_data[k + 2],
03530                                                                                                         image2_data[k + image2_nx],
03531                                                                                                         image2_data[k + 2 + image2_nx], x2, y2);
03532 
03533                         s1 += Util::square_sum(im1[l], im1[l + 1]);
03534                         s2 += Util::square_sum(im2[l], im2[l + 1]);
03535                 }
03536 
03537                 jmax = j - 1;
03538                 float sqrt_s1 = std::sqrt(s1);
03539                 float sqrt_s2 = std::sqrt(s2);
03540 
03541                 int l = 0;
03542                 for (float r = 1; r < rmax; r += 1.0f) {
03543                         int i3 = i2 + l;
03544                         im1[i3] /= sqrt_s1;
03545                         im1[i3 + 1] /= sqrt_s1;
03546                         im2[i3] /= sqrt_s2;
03547                         im2[i3 + 1] /= sqrt_s2;
03548                         l += 2;
03549                 }
03550         }
03551         float * data = get_data();
03552 
03553         switch (mode) {
03554         case 0:
03555                 for (int m1 = 0; m1 < 2; m1++) {
03556                         for (int m2 = 0; m2 < 2; m2++) {
03557 
03558                                 if (m1 == 0 && m2 == 0) {
03559                                         for (int i = 0; i < steps; i++) {
03560                                                 int i2 = i * rmax * 2;
03561                                                 for (int j = 0; j < steps; j++) {
03562                                                         int l = i + j * steps * 2;
03563                                                         int j2 = j * rmax * 2;
03564                                                         data[l] = 0;
03565                                                         for (int k = 0; k < jmax; k++) {
03566                                                                 data[l] += im1[i2 + k] * im2[j2 + k];
03567                                                         }
03568                                                 }
03569                                         }
03570                                 }
03571                                 else {
03572                                         int steps2 = steps * m2 + steps * steps * 2 * m1;
03573 
03574                                         for (int i = 0; i < steps; i++) {
03575                                                 int i2 = i * rmax * 2;
03576                                                 for (int j = 0; j < steps; j++) {
03577                                                         int j2 = j * rmax * 2;
03578                                                         int l = i + j * steps * 2 + steps2;
03579                                                         data[l] = 0;
03580 
03581                                                         for (int k = 0; k < jmax; k += 2) {
03582                                                                 i2 += k;
03583                                                                 j2 += k;
03584                                                                 data[l] += im1[i2] * im2[j2];
03585                                                                 data[l] += -im1[i2 + 1] * im2[j2 + 1];
03586                                                         }
03587                                                 }
03588                                         }
03589                                 }
03590                         }
03591                 }
03592 
03593                 break;
03594         case 1:
03595                 for (int m1 = 0; m1 < 2; m1++) {
03596                         for (int m2 = 0; m2 < 2; m2++) {
03597                                 int steps2 = steps * m2 + steps * steps * 2 * m1;
03598                                 int p1_sign = 1;
03599                                 if (m1 != m2) {
03600                                         p1_sign = -1;
03601                                 }
03602 
03603                                 for (int i = 0; i < steps; i++) {
03604                                         int i2 = i * rmax * 2;
03605 
03606                                         for (int j = 0; j < steps; j++) {
03607                                                 int j2 = j * rmax * 2;
03608 
03609                                                 int l = i + j * steps * 2 + steps2;
03610                                                 data[l] = 0;
03611                                                 float a = 0;
03612 
03613                                                 for (int k = 0; k < jmax; k += 2) {
03614                                                         i2 += k;
03615                                                         j2 += k;
03616 
03617 #ifdef  _WIN32
03618                                                         float a1 = (float) _hypot(im1[i2], im1[i2 + 1]);
03619 #else
03620                                                         float a1 = (float) hypot(im1[i2], im1[i2 + 1]);
03621 #endif  //_WIN32
03622                                                         float p1 = atan2(im1[i2 + 1], im1[i2]);
03623                                                         float p2 = atan2(im2[j2 + 1], im2[j2]);
03624 
03625                                                         data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1;
03626                                                         a += a1;
03627                                                 }
03628 
03629                                                 data[l] /= (float)(a * M_PI / 180.0f);
03630                                         }
03631                                 }
03632                         }
03633                 }
03634 
03635                 break;
03636         case 2:
03637                 for (int m1 = 0; m1 < 2; m1++) {
03638                         for (int m2 = 0; m2 < 2; m2++) {
03639                                 int steps2 = steps * m2 + steps * steps * 2 * m1;
03640 
03641                                 for (int i = 0; i < steps; i++) {
03642                                         int i2 = i * rmax * 2;
03643 
03644                                         for (int j = 0; j < steps; j++) {
03645                                                 int j2 = j * rmax * 2;
03646                                                 int l = i + j * steps * 2 + steps2;
03647                                                 data[l] = 0;
03648 
03649                                                 for (int k = 0; k < jmax; k += 2) {
03650                                                         i2 += k;
03651                                                         j2 += k;
03652 #ifdef  _WIN32
03653                                                         data[l] += (float) (_hypot(im1[i2], im1[i2 + 1]) * _hypot(im2[j2], im2[j2 + 1]));
03654 #else
03655                                                         data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1]));
03656 #endif  //_WIN32
03657                                                 }
03658                                         }
03659                                 }
03660                         }
03661                 }
03662 
03663                 break;
03664         default:
03665                 break;
03666         }
03667 
03668         if (horizontal) {
03669                 float *tmp_array = new float[ny];
03670                 for (int i = 1; i < nx; i++) {
03671                         for (int j = 0; j < ny; j++) {
03672                                 tmp_array[j] = get_value_at(i, j);
03673                         }
03674                         for (int j = 0; j < ny; j++) {
03675                                 set_value_at(i, j, 0, tmp_array[(j + i) % ny]);
03676                         }
03677                 }
03678                 if( tmp_array )
03679                 {
03680                         delete[]tmp_array;
03681                         tmp_array = 0;
03682                 }
03683         }
03684 
03685         if( im1 )
03686         {
03687                 delete[]im1;
03688                 im1 = 0;
03689         }
03690 
03691         if( im2 )
03692         {
03693                 delete im2;
03694                 im2 = 0;
03695         }
03696 
03697 
03698         image1->update();
03699         image2->update();
03700         if( image1 )
03701         {
03702                 delete image1;
03703                 image1 = 0;
03704         }
03705         if( image2 )
03706         {
03707                 delete image2;
03708                 image2 = 0;
03709         }
03710         update();
03711         EXITFUNC;
03712 }

void EMData::common_lines_real EMData image1,
EMData image2,
int  steps = 180,
bool  horizontal = false
[inherited]
 

Finds common lines between 2 real images.

Parameters:
image1 The first image.
image2 The second image.
steps 1/2 of the resolution of the map.
horizontal In horizontal way or not.
Exceptions:
NullPointerException If 'image1' or 'image2' is NULL.
ImageFormatException If 'image1' 'image2' are not same size.

Definition at line 3716 of file emdata.cpp.

References EMAN::EMData::copy(), data, EMAN::EMData::get_data(), EMAN::EMData::get_ysize(), ImageFormatException, images, NullPointerException, EMAN::EMData::set_size(), sqrt(), t, EMAN::EMData::transform(), and EMAN::EMData::update().

Referenced by main().

03718 {
03719         ENTERFUNC;
03720 
03721         if (!image1 || !image2) {
03722                 throw NullPointerException("NULL image");
03723         }
03724 
03725         if (!EMUtil::is_same_size(image1, image2)) {
03726                 throw ImageFormatException("images not same size");
03727         }
03728 
03729         int steps2 = steps * 2;
03730         int image_ny = image1->get_ysize();
03731         EMData *image1_copy = image1->copy();
03732         EMData *image2_copy = image2->copy();
03733 
03734         float *im1 = new float[steps2 * image_ny];
03735         float *im2 = new float[steps2 * image_ny];
03736 
03737         EMData *images[] = { image1_copy, image2_copy };
03738         float *ims[] = { im1, im2 };
03739 
03740         for (int m = 0; m < 2; m++) {
03741                 float *im = ims[m];
03742                 float a = M_PI / steps2;
03743                 Transform t(Dict("type","2d","alpha",-a));
03744                 for (int i = 0; i < steps2; i++) {
03745                         images[i]->transform(t);
03746                         float *data = images[i]->get_data();
03747 
03748                         for (int j = 0; j < image_ny; j++) {
03749                                 float sum = 0;
03750                                 for (int k = 0; k < image_ny; k++) {
03751                                         sum += data[j * image_ny + k];
03752                                 }
03753                                 im[i * image_ny + j] = sum;
03754                         }
03755 
03756                         float sum1 = 0;
03757                         float sum2 = 0;
03758                         for (int j = 0; j < image_ny; j++) {
03759                                 int l = i * image_ny + j;
03760                                 sum1 += im[l];
03761                                 sum2 += im[l] * im[l];
03762                         }
03763 
03764                         float mean = sum1 / image_ny;
03765                         float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1);
03766 
03767                         for (int j = 0; j < image_ny; j++) {
03768                                 int l = i * image_ny + j;
03769                                 im[l] = (im[l] - mean) / sigma;
03770                         }
03771 
03772                         images[i]->update();
03773                         a += M_PI / steps;
03774                 }
03775         }
03776 
03777         set_size(steps2, steps2, 1);
03778         float *data1 = get_data();
03779 
03780         if (horiz) {
03781                 for (int i = 0; i < steps2; i++) {
03782                         for (int j = 0, l = i; j < steps2; j++, l++) {
03783                                 if (l == steps2) {
03784                                         l = 0;
03785                                 }
03786 
03787                                 float sum = 0;
03788                                 for (int k = 0; k < image_ny; k++) {
03789                                         sum += im1[i * image_ny + k] * im2[l * image_ny + k];
03790                                 }
03791                                 data1[i + j * steps2] = sum;
03792                         }
03793                 }
03794         }
03795         else {
03796                 for (int i = 0; i < steps2; i++) {
03797                         for (int j = 0; j < steps2; j++) {
03798                                 float sum = 0;
03799                                 for (int k = 0; k < image_ny; k++) {
03800                                         sum += im1[i * image_ny + k] * im2[j * image_ny + k];
03801                                 }
03802                                 data1[i + j * steps2] = sum;
03803                         }
03804                 }
03805         }
03806 
03807         update();
03808 
03809         if( image1_copy )
03810         {
03811                 delete image1_copy;
03812                 image1_copy = 0;
03813         }
03814 
03815         if( image2_copy )
03816         {
03817                 delete image2_copy;
03818                 image2_copy = 0;
03819         }
03820 
03821         if( im1 )
03822         {
03823                 delete[]im1;
03824                 im1 = 0;
03825         }
03826 
03827         if( im2 )
03828         {
03829                 delete[]im2;
03830                 im2 = 0;
03831         }
03832         EXITFUNC;
03833 }


Generated on Fri Aug 10 16:36:39 2012 for EMAN2 by  doxygen 1.3.9.1