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 3385 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().
03387 { 03388 ENTERFUNC; 03389 03390 if (!image1 || !image2) { 03391 throw NullPointerException("NULL image"); 03392 } 03393 03394 if (mode < 0 || mode > 2) { 03395 throw OutofRangeException(0, 2, mode, "invalid mode"); 03396 } 03397 03398 if (!image1->is_complex()) { 03399 image1 = image1->do_fft(); 03400 } 03401 if (!image2->is_complex()) { 03402 image2 = image2->do_fft(); 03403 } 03404 03405 image1->ap2ri(); 03406 image2->ap2ri(); 03407 03408 if (!EMUtil::is_same_size(image1, image2)) { 03409 throw ImageFormatException("images not same sizes"); 03410 } 03411 03412 int image2_nx = image2->get_xsize(); 03413 int image2_ny = image2->get_ysize(); 03414 03415 int rmax = image2_ny / 4 - 1; 03416 int array_size = steps * rmax * 2; 03417 float *im1 = new float[array_size]; 03418 float *im2 = new float[array_size]; 03419 for (int i = 0; i < array_size; i++) { 03420 im1[i] = 0; 03421 im2[i] = 0; 03422 } 03423 03424 set_size(steps * 2, steps * 2, 1); 03425 03426 float *image1_data = image1->get_data(); 03427 float *image2_data = image2->get_data(); 03428 03429 float da = M_PI / steps; 03430 float a = -M_PI / 2.0f + da / 2.0f; 03431 int jmax = 0; 03432 03433 for (int i = 0; i < steps * 2; i += 2, a += da) { 03434 float s1 = 0; 03435 float s2 = 0; 03436 int i2 = i * rmax; 03437 int j = 0; 03438 03439 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) { 03440 float x = r * cos(a); 03441 float y = r * sin(a); 03442 03443 if (x < 0) { 03444 x = -x; 03445 y = -y; 03446 LOGERR("CCL ERROR %d, %f !\n", i, -x); 03447 } 03448 03449 int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx); 03450 int l = i2 + j; 03451 float x2 = x - floor(x); 03452 float y2 = y - floor(y); 03453 03454 im1[l] = Util::bilinear_interpolate(image1_data[k], 03455 image1_data[k + 2], 03456 image1_data[k + image2_nx], 03457 image1_data[k + 2 + image2_nx], x2, y2); 03458 03459 im2[l] = Util::bilinear_interpolate(image2_data[k], 03460 image2_data[k + 2], 03461 image2_data[k + image2_nx], 03462 image2_data[k + 2 + image2_nx], x2, y2); 03463 03464 k++; 03465 03466 im1[l + 1] = Util::bilinear_interpolate(image1_data[k], 03467 image1_data[k + 2], 03468 image1_data[k + image2_nx], 03469 image1_data[k + 2 + image2_nx], x2, y2); 03470 03471 im2[l + 1] = Util::bilinear_interpolate(image2_data[k], 03472 image2_data[k + 2], 03473 image2_data[k + image2_nx], 03474 image2_data[k + 2 + image2_nx], x2, y2); 03475 03476 s1 += Util::square_sum(im1[l], im1[l + 1]); 03477 s2 += Util::square_sum(im2[l], im2[l + 1]); 03478 } 03479 03480 jmax = j - 1; 03481 float sqrt_s1 = std::sqrt(s1); 03482 float sqrt_s2 = std::sqrt(s2); 03483 03484 int l = 0; 03485 for (float r = 1; r < rmax; r += 1.0f) { 03486 int i3 = i2 + l; 03487 im1[i3] /= sqrt_s1; 03488 im1[i3 + 1] /= sqrt_s1; 03489 im2[i3] /= sqrt_s2; 03490 im2[i3 + 1] /= sqrt_s2; 03491 l += 2; 03492 } 03493 } 03494 float * data = get_data(); 03495 03496 switch (mode) { 03497 case 0: 03498 for (int m1 = 0; m1 < 2; m1++) { 03499 for (int m2 = 0; m2 < 2; m2++) { 03500 03501 if (m1 == 0 && m2 == 0) { 03502 for (int i = 0; i < steps; i++) { 03503 int i2 = i * rmax * 2; 03504 for (int j = 0; j < steps; j++) { 03505 int l = i + j * steps * 2; 03506 int j2 = j * rmax * 2; 03507 data[l] = 0; 03508 for (int k = 0; k < jmax; k++) { 03509 data[l] += im1[i2 + k] * im2[j2 + k]; 03510 } 03511 } 03512 } 03513 } 03514 else { 03515 int steps2 = steps * m2 + steps * steps * 2 * m1; 03516 03517 for (int i = 0; i < steps; i++) { 03518 int i2 = i * rmax * 2; 03519 for (int j = 0; j < steps; j++) { 03520 int j2 = j * rmax * 2; 03521 int l = i + j * steps * 2 + steps2; 03522 data[l] = 0; 03523 03524 for (int k = 0; k < jmax; k += 2) { 03525 i2 += k; 03526 j2 += k; 03527 data[l] += im1[i2] * im2[j2]; 03528 data[l] += -im1[i2 + 1] * im2[j2 + 1]; 03529 } 03530 } 03531 } 03532 } 03533 } 03534 } 03535 03536 break; 03537 case 1: 03538 for (int m1 = 0; m1 < 2; m1++) { 03539 for (int m2 = 0; m2 < 2; m2++) { 03540 int steps2 = steps * m2 + steps * steps * 2 * m1; 03541 int p1_sign = 1; 03542 if (m1 != m2) { 03543 p1_sign = -1; 03544 } 03545 03546 for (int i = 0; i < steps; i++) { 03547 int i2 = i * rmax * 2; 03548 03549 for (int j = 0; j < steps; j++) { 03550 int j2 = j * rmax * 2; 03551 03552 int l = i + j * steps * 2 + steps2; 03553 data[l] = 0; 03554 float a = 0; 03555 03556 for (int k = 0; k < jmax; k += 2) { 03557 i2 += k; 03558 j2 += k; 03559 03560 #ifdef _WIN32 03561 float a1 = (float) _hypot(im1[i2], im1[i2 + 1]); 03562 #else 03563 float a1 = (float) hypot(im1[i2], im1[i2 + 1]); 03564 #endif //_WIN32 03565 float p1 = atan2(im1[i2 + 1], im1[i2]); 03566 float p2 = atan2(im2[j2 + 1], im2[j2]); 03567 03568 data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1; 03569 a += a1; 03570 } 03571 03572 data[l] /= (float)(a * M_PI / 180.0f); 03573 } 03574 } 03575 } 03576 } 03577 03578 break; 03579 case 2: 03580 for (int m1 = 0; m1 < 2; m1++) { 03581 for (int m2 = 0; m2 < 2; m2++) { 03582 int steps2 = steps * m2 + steps * steps * 2 * m1; 03583 03584 for (int i = 0; i < steps; i++) { 03585 int i2 = i * rmax * 2; 03586 03587 for (int j = 0; j < steps; j++) { 03588 int j2 = j * rmax * 2; 03589 int l = i + j * steps * 2 + steps2; 03590 data[l] = 0; 03591 03592 for (int k = 0; k < jmax; k += 2) { 03593 i2 += k; 03594 j2 += k; 03595 #ifdef _WIN32 03596 data[l] += (float) (_hypot(im1[i2], im1[i2 + 1]) * _hypot(im2[j2], im2[j2 + 1])); 03597 #else 03598 data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1])); 03599 #endif //_WIN32 03600 } 03601 } 03602 } 03603 } 03604 } 03605 03606 break; 03607 default: 03608 break; 03609 } 03610 03611 if (horizontal) { 03612 float *tmp_array = new float[ny]; 03613 for (int i = 1; i < nx; i++) { 03614 for (int j = 0; j < ny; j++) { 03615 tmp_array[j] = get_value_at(i, j); 03616 } 03617 for (int j = 0; j < ny; j++) { 03618 set_value_at(i, j, 0, tmp_array[(j + i) % ny]); 03619 } 03620 } 03621 if( tmp_array ) 03622 { 03623 delete[]tmp_array; 03624 tmp_array = 0; 03625 } 03626 } 03627 03628 if( im1 ) 03629 { 03630 delete[]im1; 03631 im1 = 0; 03632 } 03633 03634 if( im2 ) 03635 { 03636 delete im2; 03637 im2 = 0; 03638 } 03639 03640 03641 image1->update(); 03642 image2->update(); 03643 if( image1 ) 03644 { 03645 delete image1; 03646 image1 = 0; 03647 } 03648 if( image2 ) 03649 { 03650 delete image2; 03651 image2 = 0; 03652 } 03653 update(); 03654 EXITFUNC; 03655 }
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 3659 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().
03661 { 03662 ENTERFUNC; 03663 03664 if (!image1 || !image2) { 03665 throw NullPointerException("NULL image"); 03666 } 03667 03668 if (!EMUtil::is_same_size(image1, image2)) { 03669 throw ImageFormatException("images not same size"); 03670 } 03671 03672 int steps2 = steps * 2; 03673 int image_ny = image1->get_ysize(); 03674 EMData *image1_copy = image1->copy(); 03675 EMData *image2_copy = image2->copy(); 03676 03677 float *im1 = new float[steps2 * image_ny]; 03678 float *im2 = new float[steps2 * image_ny]; 03679 03680 EMData *images[] = { image1_copy, image2_copy }; 03681 float *ims[] = { im1, im2 }; 03682 03683 for (int m = 0; m < 2; m++) { 03684 float *im = ims[m]; 03685 float a = M_PI / steps2; 03686 Transform t(Dict("type","2d","alpha",-a)); 03687 for (int i = 0; i < steps2; i++) { 03688 images[i]->transform(t); 03689 float *data = images[i]->get_data(); 03690 03691 for (int j = 0; j < image_ny; j++) { 03692 float sum = 0; 03693 for (int k = 0; k < image_ny; k++) { 03694 sum += data[j * image_ny + k]; 03695 } 03696 im[i * image_ny + j] = sum; 03697 } 03698 03699 float sum1 = 0; 03700 float sum2 = 0; 03701 for (int j = 0; j < image_ny; j++) { 03702 int l = i * image_ny + j; 03703 sum1 += im[l]; 03704 sum2 += im[l] * im[l]; 03705 } 03706 03707 float mean = sum1 / image_ny; 03708 float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1); 03709 03710 for (int j = 0; j < image_ny; j++) { 03711 int l = i * image_ny + j; 03712 im[l] = (im[l] - mean) / sigma; 03713 } 03714 03715 images[i]->update(); 03716 a += M_PI / steps; 03717 } 03718 } 03719 03720 set_size(steps2, steps2, 1); 03721 float *data1 = get_data(); 03722 03723 if (horiz) { 03724 for (int i = 0; i < steps2; i++) { 03725 for (int j = 0, l = i; j < steps2; j++, l++) { 03726 if (l == steps2) { 03727 l = 0; 03728 } 03729 03730 float sum = 0; 03731 for (int k = 0; k < image_ny; k++) { 03732 sum += im1[i * image_ny + k] * im2[l * image_ny + k]; 03733 } 03734 data1[i + j * steps2] = sum; 03735 } 03736 } 03737 } 03738 else { 03739 for (int i = 0; i < steps2; i++) { 03740 for (int j = 0; j < steps2; j++) { 03741 float sum = 0; 03742 for (int k = 0; k < image_ny; k++) { 03743 sum += im1[i * image_ny + k] * im2[j * image_ny + k]; 03744 } 03745 data1[i + j * steps2] = sum; 03746 } 03747 } 03748 } 03749 03750 update(); 03751 03752 if( image1_copy ) 03753 { 03754 delete image1_copy; 03755 image1_copy = 0; 03756 } 03757 03758 if( image2_copy ) 03759 { 03760 delete image2_copy; 03761 image2_copy = 0; 03762 } 03763 03764 if( im1 ) 03765 { 03766 delete[]im1; 03767 im1 = 0; 03768 } 03769 03770 if( im2 ) 03771 { 03772 delete[]im2; 03773 im2 = 0; 03774 } 03775 EXITFUNC; 03776 }