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

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

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

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


Generated on Thu May 3 10:08:19 2012 for EMAN2 by  doxygen 1.4.7