EMAN2
Functions
code complete but contains bugs

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. More...
 
void EMAN::EMData::common_lines_real (EMData *image1, EMData *image2, int steps=180, bool horizontal=false)
 Finds common lines between 2 real images. More...
 

Detailed Description

Function Documentation

◆ common_lines()

void EMData::common_lines ( EMData *  image1,
EMData *  image2,
int  mode = 0,
int  steps = 180,
bool  horizontal = false 
)

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.

Parameters
image1The first complex image.
image2The second complex image.
modeEither 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.
steps1/2 of the resolution of the map.
horizontalIn horizontal way or not.
Exceptions
NullPointerExceptionIf 'image1' or 'image2' is NULL.
OutofRangeExceptionIf 'mode' is invalid.
ImageFormatExceptionIf 'image1' 'image2' are not same size.

Definition at line 3632 of file emdata.cpp.

3634{
3635 ENTERFUNC;
3636
3637 if (!image1 || !image2) {
3638 throw NullPointerException("NULL image");
3639 }
3640
3641 if (mode < 0 || mode > 2) {
3642 throw OutofRangeException(0, 2, mode, "invalid mode");
3643 }
3644
3645 if (!image1->is_complex()) {
3646 image1 = image1->do_fft();
3647 }
3648 if (!image2->is_complex()) {
3649 image2 = image2->do_fft();
3650 }
3651
3652 image1->ap2ri();
3653 image2->ap2ri();
3654
3655 if (!EMUtil::is_same_size(image1, image2)) {
3656 throw ImageFormatException("images not same sizes");
3657 }
3658
3659 int image2_nx = image2->get_xsize();
3660 int image2_ny = image2->get_ysize();
3661
3662 int rmax = image2_ny / 4 - 1;
3663 int array_size = steps * rmax * 2;
3664 float *im1 = new float[array_size];
3665 float *im2 = new float[array_size];
3666 for (int i = 0; i < array_size; i++) {
3667 im1[i] = 0;
3668 im2[i] = 0;
3669 }
3670
3671 set_size(steps * 2, steps * 2, 1);
3672
3673 float *image1_data = image1->get_data();
3674 float *image2_data = image2->get_data();
3675
3676 float da = M_PI / steps;
3677 float a = -M_PI / 2.0f + da / 2.0f;
3678 int jmax = 0;
3679
3680 for (int i = 0; i < steps * 2; i += 2, a += da) {
3681 float s1 = 0;
3682 float s2 = 0;
3683 int i2 = i * rmax;
3684 int j = 0;
3685
3686 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) {
3687 float x = r * cos(a);
3688 float y = r * sin(a);
3689
3690 if (x < 0) {
3691 x = -x;
3692 y = -y;
3693 LOGERR("CCL ERROR %d, %f !\n", i, -x);
3694 }
3695
3696 int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx);
3697 int l = i2 + j;
3698 float x2 = x - floor(x);
3699 float y2 = y - floor(y);
3700
3701 im1[l] = Util::bilinear_interpolate(image1_data[k],
3702 image1_data[k + 2],
3703 image1_data[k + image2_nx],
3704 image1_data[k + 2 + image2_nx], x2, y2);
3705
3706 im2[l] = Util::bilinear_interpolate(image2_data[k],
3707 image2_data[k + 2],
3708 image2_data[k + image2_nx],
3709 image2_data[k + 2 + image2_nx], x2, y2);
3710
3711 k++;
3712
3713 im1[l + 1] = Util::bilinear_interpolate(image1_data[k],
3714 image1_data[k + 2],
3715 image1_data[k + image2_nx],
3716 image1_data[k + 2 + image2_nx], x2, y2);
3717
3718 im2[l + 1] = Util::bilinear_interpolate(image2_data[k],
3719 image2_data[k + 2],
3720 image2_data[k + image2_nx],
3721 image2_data[k + 2 + image2_nx], x2, y2);
3722
3723 s1 += Util::square_sum(im1[l], im1[l + 1]);
3724 s2 += Util::square_sum(im2[l], im2[l + 1]);
3725 }
3726
3727 jmax = j - 1;
3728 float sqrt_s1 = std::sqrt(s1);
3729 float sqrt_s2 = std::sqrt(s2);
3730
3731 int l = 0;
3732 for (float r = 1; r < rmax; r += 1.0f) {
3733 int i3 = i2 + l;
3734 im1[i3] /= sqrt_s1;
3735 im1[i3 + 1] /= sqrt_s1;
3736 im2[i3] /= sqrt_s2;
3737 im2[i3 + 1] /= sqrt_s2;
3738 l += 2;
3739 }
3740 }
3741 float * data = get_data();
3742
3743 switch (mode) {
3744 case 0:
3745 for (int m1 = 0; m1 < 2; m1++) {
3746 for (int m2 = 0; m2 < 2; m2++) {
3747
3748 if (m1 == 0 && m2 == 0) {
3749 for (int i = 0; i < steps; i++) {
3750 int i2 = i * rmax * 2;
3751 for (int j = 0; j < steps; j++) {
3752 int l = i + j * steps * 2;
3753 int j2 = j * rmax * 2;
3754 data[l] = 0;
3755 for (int k = 0; k < jmax; k++) {
3756 data[l] += im1[i2 + k] * im2[j2 + k];
3757 }
3758 }
3759 }
3760 }
3761 else {
3762 int steps2 = steps * m2 + steps * steps * 2 * m1;
3763
3764 for (int i = 0; i < steps; i++) {
3765 int i2 = i * rmax * 2;
3766 for (int j = 0; j < steps; j++) {
3767 int j2 = j * rmax * 2;
3768 int l = i + j * steps * 2 + steps2;
3769 data[l] = 0;
3770
3771 for (int k = 0; k < jmax; k += 2) {
3772 i2 += k;
3773 j2 += k;
3774 data[l] += im1[i2] * im2[j2];
3775 data[l] += -im1[i2 + 1] * im2[j2 + 1];
3776 }
3777 }
3778 }
3779 }
3780 }
3781 }
3782
3783 break;
3784 case 1:
3785 for (int m1 = 0; m1 < 2; m1++) {
3786 for (int m2 = 0; m2 < 2; m2++) {
3787 int steps2 = steps * m2 + steps * steps * 2 * m1;
3788 int p1_sign = 1;
3789 if (m1 != m2) {
3790 p1_sign = -1;
3791 }
3792
3793 for (int i = 0; i < steps; i++) {
3794 int i2 = i * rmax * 2;
3795
3796 for (int j = 0; j < steps; j++) {
3797 int j2 = j * rmax * 2;
3798
3799 int l = i + j * steps * 2 + steps2;
3800 data[l] = 0;
3801 float a = 0;
3802
3803 for (int k = 0; k < jmax; k += 2) {
3804 i2 += k;
3805 j2 += k;
3806
3807 float a1 = (float) hypot(im1[i2], im1[i2 + 1]);
3808 float p1 = atan2(im1[i2 + 1], im1[i2]);
3809 float p2 = atan2(im2[j2 + 1], im2[j2]);
3810
3811 data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1;
3812 a += a1;
3813 }
3814
3815 data[l] /= (float)(a * M_PI / 180.0f);
3816 }
3817 }
3818 }
3819 }
3820
3821 break;
3822 case 2:
3823 for (int m1 = 0; m1 < 2; m1++) {
3824 for (int m2 = 0; m2 < 2; m2++) {
3825 int steps2 = steps * m2 + steps * steps * 2 * m1;
3826
3827 for (int i = 0; i < steps; i++) {
3828 int i2 = i * rmax * 2;
3829
3830 for (int j = 0; j < steps; j++) {
3831 int j2 = j * rmax * 2;
3832 int l = i + j * steps * 2 + steps2;
3833 data[l] = 0;
3834
3835 for (int k = 0; k < jmax; k += 2) {
3836 i2 += k;
3837 j2 += k;
3838 data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1]));
3839 }
3840 }
3841 }
3842 }
3843 }
3844
3845 break;
3846 default:
3847 break;
3848 }
3849
3850 if (horizontal) {
3851 float *tmp_array = new float[ny];
3852 for (int i = 1; i < nx; i++) {
3853 for (int j = 0; j < ny; j++) {
3854 tmp_array[j] = get_value_at(i, j);
3855 }
3856 for (int j = 0; j < ny; j++) {
3857 set_value_at(i, j, 0, tmp_array[(j + i) % ny]);
3858 }
3859 }
3860 if( tmp_array )
3861 {
3862 delete[]tmp_array;
3863 tmp_array = 0;
3864 }
3865 }
3866
3867 if( im1 )
3868 {
3869 delete[]im1;
3870 im1 = 0;
3871 }
3872
3873 if( im2 )
3874 {
3875 delete[] im2;
3876 im2 = 0;
3877 }
3878
3879
3880 image1->update();
3881 image2->update();
3882 if( image1 )
3883 {
3884 delete image1;
3885 image1 = 0;
3886 }
3887 if( image2 )
3888 {
3889 delete image2;
3890 image2 = 0;
3891 }
3892 update();
3893 EXITFUNC;
3894}
EMAN::EMData::nx
int nx
image size
Definition: emdata.h:848
EMAN::EMUtil::is_same_size
static bool is_same_size(const EMData *image1, const EMData *image2)
Check whether two EMData images are of the same size.
Definition: emutil.cpp:1224
EMAN::Util::angle_sub_2pi
static float angle_sub_2pi(float x, float y)
Calculate the difference of 2 angles and makes the equivalent result to be less than Pi.
Definition: util.h:1066
EMAN::Util::bilinear_interpolate
static float bilinear_interpolate(float p1, float p2, float p3, float p4, float t, float u)
Calculate bilinear interpolation.
Definition: util.h:543
EMAN::Util::square_sum
static float square_sum(float x, float y)
Calcuate (x*x + y*y).
Definition: util.h:764
set_value_at
void set_value_at(Vec3i loc, float val)
set_value_at with Vec3i
Definition: emdata_core.h:552
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
get_value_at
float get_value_at(int x, int y, int z) const
Get the pixel density value at coordinates (x,y,z).
Definition: emdata_core.h:221
set_size
void set_size(int nx, int ny=1, int nz=1, bool noalloc=false)
Resize this EMData's main board memory pointer.
update
void update()
Mark EMData as changed, statistics, etc will be updated at need.
Definition: emdata_metadata.h:138
get_data
float * get_data() const
Get the image pixel density data in a 1D float array.
Definition: emdata_metadata.h:78
ImageFormatException
#define ImageFormatException(desc)
Definition: exception.h:147
OutofRangeException
#define OutofRangeException(low, high, input, objname)
Definition: exception.h:334
NullPointerException
#define NullPointerException(desc)
Definition: exception.h:241
LOGERR
#define LOGERR
Definition: log.h:51
ENTERFUNC
#define ENTERFUNC
Definition: log.h:48
EXITFUNC
#define EXITFUNC
Definition: log.h:49
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517

References EMAN::Util::angle_sub_2pi(), EMAN::Util::bilinear_interpolate(), ENTERFUNC, EXITFUNC, get_data(), get_value_at(), ImageFormatException, EMAN::EMUtil::is_same_size(), LOGERR, NullPointerException, EMAN::EMData::nx, EMAN::EMData::ny, OutofRangeException, set_size(), set_value_at(), sqrt(), EMAN::Util::square_sum(), update(), x, and y.

◆ common_lines_real()

void EMData::common_lines_real ( EMData *  image1,
EMData *  image2,
int  steps = 180,
bool  horizontal = false 
)

Finds common lines between 2 real images.

Parameters
image1The first image.
image2The second image.
steps1/2 of the resolution of the map.
horizontalIn horizontal way or not.
Exceptions
NullPointerExceptionIf 'image1' or 'image2' is NULL.
ImageFormatExceptionIf 'image1' 'image2' are not same size.

Definition at line 3898 of file emdata.cpp.

3900{
3901 ENTERFUNC;
3902
3903 if (!image1 || !image2) {
3904 throw NullPointerException("NULL image");
3905 }
3906
3907 if (!EMUtil::is_same_size(image1, image2)) {
3908 throw ImageFormatException("images not same size");
3909 }
3910
3911 int steps2 = steps * 2;
3912 int image_ny = image1->get_ysize();
3913 EMData *image1_copy = image1->copy();
3914 EMData *image2_copy = image2->copy();
3915
3916 float *im1 = new float[steps2 * image_ny];
3917 float *im2 = new float[steps2 * image_ny];
3918
3919 EMData *images[] = { image1_copy, image2_copy };
3920 float *ims[] = { im1, im2 };
3921
3922 for (int m = 0; m < 2; m++) {
3923 float *im = ims[m];
3924 float a = M_PI / steps2;
3925 Transform t(Dict("type","2d","alpha",-a));
3926 for (int i = 0; i < steps2; i++) {
3927 images[i]->transform(t);
3928 float *data = images[i]->get_data();
3929
3930 for (int j = 0; j < image_ny; j++) {
3931 float sum = 0;
3932 for (int k = 0; k < image_ny; k++) {
3933 sum += data[j * image_ny + k];
3934 }
3935 im[i * image_ny + j] = sum;
3936 }
3937
3938 float sum1 = 0;
3939 float sum2 = 0;
3940 for (int j = 0; j < image_ny; j++) {
3941 int l = i * image_ny + j;
3942 sum1 += im[l];
3943 sum2 += im[l] * im[l];
3944 }
3945
3946 float mean = sum1 / image_ny;
3947 float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1);
3948
3949 for (int j = 0; j < image_ny; j++) {
3950 int l = i * image_ny + j;
3951 im[l] = (im[l] - mean) / sigma;
3952 }
3953
3954 images[i]->update();
3955 a += M_PI / steps;
3956 }
3957 }
3958
3959 set_size(steps2, steps2, 1);
3960 float *data1 = get_data();
3961
3962 if (horiz) {
3963 for (int i = 0; i < steps2; i++) {
3964 for (int j = 0, l = i; j < steps2; j++, l++) {
3965 if (l == steps2) {
3966 l = 0;
3967 }
3968
3969 float sum = 0;
3970 for (int k = 0; k < image_ny; k++) {
3971 sum += im1[i * image_ny + k] * im2[l * image_ny + k];
3972 }
3973 data1[i + j * steps2] = sum;
3974 }
3975 }
3976 }
3977 else {
3978 for (int i = 0; i < steps2; i++) {
3979 for (int j = 0; j < steps2; j++) {
3980 float sum = 0;
3981 for (int k = 0; k < image_ny; k++) {
3982 sum += im1[i * image_ny + k] * im2[j * image_ny + k];
3983 }
3984 data1[i + j * steps2] = sum;
3985 }
3986 }
3987 }
3988
3989 update();
3990
3991 if( image1_copy )
3992 {
3993 delete image1_copy;
3994 image1_copy = 0;
3995 }
3996
3997 if( image2_copy )
3998 {
3999 delete image2_copy;
4000 image2_copy = 0;
4001 }
4002
4003 if( im1 )
4004 {
4005 delete[]im1;
4006 im1 = 0;
4007 }
4008
4009 if( im2 )
4010 {
4011 delete[]im2;
4012 im2 = 0;
4013 }
4014 EXITFUNC;
4015}
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
images
#define images(i, j, k)
Definition: projector.cpp:1897

References ENTERFUNC, EXITFUNC, get_data(), ImageFormatException, images, EMAN::EMUtil::is_same_size(), NullPointerException, set_size(), sqrt(), and update().

Generated by doxygen 1.9.2