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 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.
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 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 }