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 3460 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(). 03462 { 03463 ENTERFUNC; 03464 03465 if (!image1 || !image2) { 03466 throw NullPointerException("NULL image"); 03467 } 03468 03469 if (mode < 0 || mode > 2) { 03470 throw OutofRangeException(0, 2, mode, "invalid mode"); 03471 } 03472 03473 if (!image1->is_complex()) { 03474 image1 = image1->do_fft(); 03475 } 03476 if (!image2->is_complex()) { 03477 image2 = image2->do_fft(); 03478 } 03479 03480 image1->ap2ri(); 03481 image2->ap2ri(); 03482 03483 if (!EMUtil::is_same_size(image1, image2)) { 03484 throw ImageFormatException("images not same sizes"); 03485 } 03486 03487 int image2_nx = image2->get_xsize(); 03488 int image2_ny = image2->get_ysize(); 03489 03490 int rmax = image2_ny / 4 - 1; 03491 int array_size = steps * rmax * 2; 03492 float *im1 = new float[array_size]; 03493 float *im2 = new float[array_size]; 03494 for (int i = 0; i < array_size; i++) { 03495 im1[i] = 0; 03496 im2[i] = 0; 03497 } 03498 03499 set_size(steps * 2, steps * 2, 1); 03500 03501 float *image1_data = image1->get_data(); 03502 float *image2_data = image2->get_data(); 03503 03504 float da = M_PI / steps; 03505 float a = -M_PI / 2.0f + da / 2.0f; 03506 int jmax = 0; 03507 03508 for (int i = 0; i < steps * 2; i += 2, a += da) { 03509 float s1 = 0; 03510 float s2 = 0; 03511 int i2 = i * rmax; 03512 int j = 0; 03513 03514 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) { 03515 float x = r * cos(a); 03516 float y = r * sin(a); 03517 03518 if (x < 0) { 03519 x = -x; 03520 y = -y; 03521 LOGERR("CCL ERROR %d, %f !\n", i, -x); 03522 } 03523 03524 int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx); 03525 int l = i2 + j; 03526 float x2 = x - floor(x); 03527 float y2 = y - floor(y); 03528 03529 im1[l] = Util::bilinear_interpolate(image1_data[k], 03530 image1_data[k + 2], 03531 image1_data[k + image2_nx], 03532 image1_data[k + 2 + image2_nx], x2, y2); 03533 03534 im2[l] = Util::bilinear_interpolate(image2_data[k], 03535 image2_data[k + 2], 03536 image2_data[k + image2_nx], 03537 image2_data[k + 2 + image2_nx], x2, y2); 03538 03539 k++; 03540 03541 im1[l + 1] = Util::bilinear_interpolate(image1_data[k], 03542 image1_data[k + 2], 03543 image1_data[k + image2_nx], 03544 image1_data[k + 2 + image2_nx], x2, y2); 03545 03546 im2[l + 1] = Util::bilinear_interpolate(image2_data[k], 03547 image2_data[k + 2], 03548 image2_data[k + image2_nx], 03549 image2_data[k + 2 + image2_nx], x2, y2); 03550 03551 s1 += Util::square_sum(im1[l], im1[l + 1]); 03552 s2 += Util::square_sum(im2[l], im2[l + 1]); 03553 } 03554 03555 jmax = j - 1; 03556 float sqrt_s1 = std::sqrt(s1); 03557 float sqrt_s2 = std::sqrt(s2); 03558 03559 int l = 0; 03560 for (float r = 1; r < rmax; r += 1.0f) { 03561 int i3 = i2 + l; 03562 im1[i3] /= sqrt_s1; 03563 im1[i3 + 1] /= sqrt_s1; 03564 im2[i3] /= sqrt_s2; 03565 im2[i3 + 1] /= sqrt_s2; 03566 l += 2; 03567 } 03568 } 03569 float * data = get_data(); 03570 03571 switch (mode) { 03572 case 0: 03573 for (int m1 = 0; m1 < 2; m1++) { 03574 for (int m2 = 0; m2 < 2; m2++) { 03575 03576 if (m1 == 0 && m2 == 0) { 03577 for (int i = 0; i < steps; i++) { 03578 int i2 = i * rmax * 2; 03579 for (int j = 0; j < steps; j++) { 03580 int l = i + j * steps * 2; 03581 int j2 = j * rmax * 2; 03582 data[l] = 0; 03583 for (int k = 0; k < jmax; k++) { 03584 data[l] += im1[i2 + k] * im2[j2 + k]; 03585 } 03586 } 03587 } 03588 } 03589 else { 03590 int steps2 = steps * m2 + steps * steps * 2 * m1; 03591 03592 for (int i = 0; i < steps; i++) { 03593 int i2 = i * rmax * 2; 03594 for (int j = 0; j < steps; j++) { 03595 int j2 = j * rmax * 2; 03596 int l = i + j * steps * 2 + steps2; 03597 data[l] = 0; 03598 03599 for (int k = 0; k < jmax; k += 2) { 03600 i2 += k; 03601 j2 += k; 03602 data[l] += im1[i2] * im2[j2]; 03603 data[l] += -im1[i2 + 1] * im2[j2 + 1]; 03604 } 03605 } 03606 } 03607 } 03608 } 03609 } 03610 03611 break; 03612 case 1: 03613 for (int m1 = 0; m1 < 2; m1++) { 03614 for (int m2 = 0; m2 < 2; m2++) { 03615 int steps2 = steps * m2 + steps * steps * 2 * m1; 03616 int p1_sign = 1; 03617 if (m1 != m2) { 03618 p1_sign = -1; 03619 } 03620 03621 for (int i = 0; i < steps; i++) { 03622 int i2 = i * rmax * 2; 03623 03624 for (int j = 0; j < steps; j++) { 03625 int j2 = j * rmax * 2; 03626 03627 int l = i + j * steps * 2 + steps2; 03628 data[l] = 0; 03629 float a = 0; 03630 03631 for (int k = 0; k < jmax; k += 2) { 03632 i2 += k; 03633 j2 += k; 03634 03635 #ifdef _WIN32 03636 float a1 = (float) _hypot(im1[i2], im1[i2 + 1]); 03637 #else 03638 float a1 = (float) hypot(im1[i2], im1[i2 + 1]); 03639 #endif //_WIN32 03640 float p1 = atan2(im1[i2 + 1], im1[i2]); 03641 float p2 = atan2(im2[j2 + 1], im2[j2]); 03642 03643 data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1; 03644 a += a1; 03645 } 03646 03647 data[l] /= (float)(a * M_PI / 180.0f); 03648 } 03649 } 03650 } 03651 } 03652 03653 break; 03654 case 2: 03655 for (int m1 = 0; m1 < 2; m1++) { 03656 for (int m2 = 0; m2 < 2; m2++) { 03657 int steps2 = steps * m2 + steps * steps * 2 * m1; 03658 03659 for (int i = 0; i < steps; i++) { 03660 int i2 = i * rmax * 2; 03661 03662 for (int j = 0; j < steps; j++) { 03663 int j2 = j * rmax * 2; 03664 int l = i + j * steps * 2 + steps2; 03665 data[l] = 0; 03666 03667 for (int k = 0; k < jmax; k += 2) { 03668 i2 += k; 03669 j2 += k; 03670 #ifdef _WIN32 03671 data[l] += (float) (_hypot(im1[i2], im1[i2 + 1]) * _hypot(im2[j2], im2[j2 + 1])); 03672 #else 03673 data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1])); 03674 #endif //_WIN32 03675 } 03676 } 03677 } 03678 } 03679 } 03680 03681 break; 03682 default: 03683 break; 03684 } 03685 03686 if (horizontal) { 03687 float *tmp_array = new float[ny]; 03688 for (int i = 1; i < nx; i++) { 03689 for (int j = 0; j < ny; j++) { 03690 tmp_array[j] = get_value_at(i, j); 03691 } 03692 for (int j = 0; j < ny; j++) { 03693 set_value_at(i, j, 0, tmp_array[(j + i) % ny]); 03694 } 03695 } 03696 if( tmp_array ) 03697 { 03698 delete[]tmp_array; 03699 tmp_array = 0; 03700 } 03701 } 03702 03703 if( im1 ) 03704 { 03705 delete[]im1; 03706 im1 = 0; 03707 } 03708 03709 if( im2 ) 03710 { 03711 delete im2; 03712 im2 = 0; 03713 } 03714 03715 03716 image1->update(); 03717 image2->update(); 03718 if( image1 ) 03719 { 03720 delete image1; 03721 image1 = 0; 03722 } 03723 if( image2 ) 03724 { 03725 delete image2; 03726 image2 = 0; 03727 } 03728 update(); 03729 EXITFUNC; 03730 }
|
|
Finds common lines between 2 real images.
Definition at line 3734 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(). 03736 { 03737 ENTERFUNC; 03738 03739 if (!image1 || !image2) { 03740 throw NullPointerException("NULL image"); 03741 } 03742 03743 if (!EMUtil::is_same_size(image1, image2)) { 03744 throw ImageFormatException("images not same size"); 03745 } 03746 03747 int steps2 = steps * 2; 03748 int image_ny = image1->get_ysize(); 03749 EMData *image1_copy = image1->copy(); 03750 EMData *image2_copy = image2->copy(); 03751 03752 float *im1 = new float[steps2 * image_ny]; 03753 float *im2 = new float[steps2 * image_ny]; 03754 03755 EMData *images[] = { image1_copy, image2_copy }; 03756 float *ims[] = { im1, im2 }; 03757 03758 for (int m = 0; m < 2; m++) { 03759 float *im = ims[m]; 03760 float a = M_PI / steps2; 03761 Transform t(Dict("type","2d","alpha",-a)); 03762 for (int i = 0; i < steps2; i++) { 03763 images[i]->transform(t); 03764 float *data = images[i]->get_data(); 03765 03766 for (int j = 0; j < image_ny; j++) { 03767 float sum = 0; 03768 for (int k = 0; k < image_ny; k++) { 03769 sum += data[j * image_ny + k]; 03770 } 03771 im[i * image_ny + j] = sum; 03772 } 03773 03774 float sum1 = 0; 03775 float sum2 = 0; 03776 for (int j = 0; j < image_ny; j++) { 03777 int l = i * image_ny + j; 03778 sum1 += im[l]; 03779 sum2 += im[l] * im[l]; 03780 } 03781 03782 float mean = sum1 / image_ny; 03783 float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1); 03784 03785 for (int j = 0; j < image_ny; j++) { 03786 int l = i * image_ny + j; 03787 im[l] = (im[l] - mean) / sigma; 03788 } 03789 03790 images[i]->update(); 03791 a += M_PI / steps; 03792 } 03793 } 03794 03795 set_size(steps2, steps2, 1); 03796 float *data1 = get_data(); 03797 03798 if (horiz) { 03799 for (int i = 0; i < steps2; i++) { 03800 for (int j = 0, l = i; j < steps2; j++, l++) { 03801 if (l == steps2) { 03802 l = 0; 03803 } 03804 03805 float sum = 0; 03806 for (int k = 0; k < image_ny; k++) { 03807 sum += im1[i * image_ny + k] * im2[l * image_ny + k]; 03808 } 03809 data1[i + j * steps2] = sum; 03810 } 03811 } 03812 } 03813 else { 03814 for (int i = 0; i < steps2; i++) { 03815 for (int j = 0; j < steps2; j++) { 03816 float sum = 0; 03817 for (int k = 0; k < image_ny; k++) { 03818 sum += im1[i * image_ny + k] * im2[j * image_ny + k]; 03819 } 03820 data1[i + j * steps2] = sum; 03821 } 03822 } 03823 } 03824 03825 update(); 03826 03827 if( image1_copy ) 03828 { 03829 delete image1_copy; 03830 image1_copy = 0; 03831 } 03832 03833 if( image2_copy ) 03834 { 03835 delete image2_copy; 03836 image2_copy = 0; 03837 } 03838 03839 if( im1 ) 03840 { 03841 delete[]im1; 03842 im1 = 0; 03843 } 03844 03845 if( im2 ) 03846 { 03847 delete[]im2; 03848 im2 = 0; 03849 } 03850 EXITFUNC; 03851 }
|