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