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. | |
|
||||||||||||||||||||||||
|
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.
Definition at line 3415 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(). 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 }
|
|
||||||||||||||||||||
|
Finds common lines between 2 real images.
Definition at line 3689 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(). 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 }
|
1.3.9.1