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 3460 of file emdata.cpp.

References EMAN::Util::angle_sub_2pi(), EMAN::EMData::ap2ri(), EMAN::Util::bilinear_interpolate(), data, EMAN::EMData::do_fft(), ENTERFUNC, EXITFUNC, EMAN::EMData::get_data(), EMAN::EMData::get_value_at(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), ImageFormatException, EMAN::EMData::is_complex(), EMAN::EMUtil::is_same_size(), LOGERR, NullPointerException, EMAN::EMData::nx, EMAN::EMData::ny, OutofRangeException, EMAN::EMData::set_size(), EMAN::EMData::set_value_at(), sqrt(), EMAN::Util::square_sum(), EMAN::EMData::update(), x, and y.

Referenced by main().

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

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 3734 of file emdata.cpp.

References EMAN::EMData::copy(), data, ENTERFUNC, EXITFUNC, EMAN::EMData::get_data(), ImageFormatException, images, EMAN::EMUtil::is_same_size(), mean(), NullPointerException, EMAN::EMData::set_size(), sqrt(), t, EMAN::EMData::transform(), and EMAN::EMData::update().

Referenced by main().

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


Generated on Tue Jun 11 12:42:12 2013 for EMAN2 by  doxygen 1.4.7