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 3415 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().

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

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 3689 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().

03691 {
03692         ENTERFUNC;
03693 
03694         if (!image1 || !image2) {
03695                 throw NullPointerException("NULL image");
03696         }
03697 
03698         if (!EMUtil::is_same_size(image1, image2)) {
03699                 throw ImageFormatException("images not same size");
03700         }
03701 
03702         int steps2 = steps * 2;
03703         int image_ny = image1->get_ysize();
03704         EMData *image1_copy = image1->copy();
03705         EMData *image2_copy = image2->copy();
03706 
03707         float *im1 = new float[steps2 * image_ny];
03708         float *im2 = new float[steps2 * image_ny];
03709 
03710         EMData *images[] = { image1_copy, image2_copy };
03711         float *ims[] = { im1, im2 };
03712 
03713         for (int m = 0; m < 2; m++) {
03714                 float *im = ims[m];
03715                 float a = M_PI / steps2;
03716                 Transform t(Dict("type","2d","alpha",-a));
03717                 for (int i = 0; i < steps2; i++) {
03718                         images[i]->transform(t);
03719                         float *data = images[i]->get_data();
03720 
03721                         for (int j = 0; j < image_ny; j++) {
03722                                 float sum = 0;
03723                                 for (int k = 0; k < image_ny; k++) {
03724                                         sum += data[j * image_ny + k];
03725                                 }
03726                                 im[i * image_ny + j] = sum;
03727                         }
03728 
03729                         float sum1 = 0;
03730                         float sum2 = 0;
03731                         for (int j = 0; j < image_ny; j++) {
03732                                 int l = i * image_ny + j;
03733                                 sum1 += im[l];
03734                                 sum2 += im[l] * im[l];
03735                         }
03736 
03737                         float mean = sum1 / image_ny;
03738                         float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1);
03739 
03740                         for (int j = 0; j < image_ny; j++) {
03741                                 int l = i * image_ny + j;
03742                                 im[l] = (im[l] - mean) / sigma;
03743                         }
03744 
03745                         images[i]->update();
03746                         a += M_PI / steps;
03747                 }
03748         }
03749 
03750         set_size(steps2, steps2, 1);
03751         float *data1 = get_data();
03752 
03753         if (horiz) {
03754                 for (int i = 0; i < steps2; i++) {
03755                         for (int j = 0, l = i; j < steps2; j++, l++) {
03756                                 if (l == steps2) {
03757                                         l = 0;
03758                                 }
03759 
03760                                 float sum = 0;
03761                                 for (int k = 0; k < image_ny; k++) {
03762                                         sum += im1[i * image_ny + k] * im2[l * image_ny + k];
03763                                 }
03764                                 data1[i + j * steps2] = sum;
03765                         }
03766                 }
03767         }
03768         else {
03769                 for (int i = 0; i < steps2; i++) {
03770                         for (int j = 0; j < steps2; j++) {
03771                                 float sum = 0;
03772                                 for (int k = 0; k < image_ny; k++) {
03773                                         sum += im1[i * image_ny + k] * im2[j * image_ny + k];
03774                                 }
03775                                 data1[i + j * steps2] = sum;
03776                         }
03777                 }
03778         }
03779 
03780         update();
03781 
03782         if( image1_copy )
03783         {
03784                 delete image1_copy;
03785                 image1_copy = 0;
03786         }
03787 
03788         if( image2_copy )
03789         {
03790                 delete image2_copy;
03791                 image2_copy = 0;
03792         }
03793 
03794         if( im1 )
03795         {
03796                 delete[]im1;
03797                 im1 = 0;
03798         }
03799 
03800         if( im2 )
03801         {
03802                 delete[]im2;
03803                 im2 = 0;
03804         }
03805         EXITFUNC;
03806 }


Generated on Tue Jul 12 13:47:31 2011 for EMAN2 by  doxygen 1.4.7