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. |
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.
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. |
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::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().
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.
image1 | The first image. | |
image2 | The second image. | |
steps | 1/2 of the resolution of the map. | |
horizontal | In horizontal way or not. |
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, 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().
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 }