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

Function Documentation

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 3479 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(), EMAN::EMData::sqrt(), EMAN::Util::square_sum(), EMAN::EMData::update(), x, and y.

Referenced by main().

{
        ENTERFUNC;

        if (!image1 || !image2) {
                throw NullPointerException("NULL image");
        }

        if (mode < 0 || mode > 2) {
                throw OutofRangeException(0, 2, mode, "invalid mode");
        }

        if (!image1->is_complex()) {
                image1 = image1->do_fft();
        }
        if (!image2->is_complex()) {
                image2 = image2->do_fft();
        }

        image1->ap2ri();
        image2->ap2ri();

        if (!EMUtil::is_same_size(image1, image2)) {
                throw ImageFormatException("images not same sizes");
        }

        int image2_nx = image2->get_xsize();
        int image2_ny = image2->get_ysize();

        int rmax = image2_ny / 4 - 1;
        int array_size = steps * rmax * 2;
        float *im1 = new float[array_size];
        float *im2 = new float[array_size];
        for (int i = 0; i < array_size; i++) {
                im1[i] = 0;
                im2[i] = 0;
        }

        set_size(steps * 2, steps * 2, 1);

        float *image1_data = image1->get_data();
        float *image2_data = image2->get_data();

        float da = M_PI / steps;
        float a = -M_PI / 2.0f + da / 2.0f;
        int jmax = 0;

        for (int i = 0; i < steps * 2; i += 2, a += da) {
                float s1 = 0;
                float s2 = 0;
                int i2 = i * rmax;
                int j = 0;

                for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) {
                        float x = r * cos(a);
                        float y = r * sin(a);

                        if (x < 0) {
                                x = -x;
                                y = -y;
                                LOGERR("CCL ERROR %d, %f !\n", i, -x);
                        }

                        int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx);
                        int l = i2 + j;
                        float x2 = x - floor(x);
                        float y2 = y - floor(y);

                        im1[l] = Util::bilinear_interpolate(image1_data[k],
                                                                                                image1_data[k + 2],
                                                                                                image1_data[k + image2_nx],
                                                                                                image1_data[k + 2 + image2_nx], x2, y2);

                        im2[l] = Util::bilinear_interpolate(image2_data[k],
                                                                                                image2_data[k + 2],
                                                                                                image2_data[k + image2_nx],
                                                                                                image2_data[k + 2 + image2_nx], x2, y2);

                        k++;

                        im1[l + 1] = Util::bilinear_interpolate(image1_data[k],
                                                                                                        image1_data[k + 2],
                                                                                                        image1_data[k + image2_nx],
                                                                                                        image1_data[k + 2 + image2_nx], x2, y2);

                        im2[l + 1] = Util::bilinear_interpolate(image2_data[k],
                                                                                                        image2_data[k + 2],
                                                                                                        image2_data[k + image2_nx],
                                                                                                        image2_data[k + 2 + image2_nx], x2, y2);

                        s1 += Util::square_sum(im1[l], im1[l + 1]);
                        s2 += Util::square_sum(im2[l], im2[l + 1]);
                }

                jmax = j - 1;
                float sqrt_s1 = std::sqrt(s1);
                float sqrt_s2 = std::sqrt(s2);

                int l = 0;
                for (float r = 1; r < rmax; r += 1.0f) {
                        int i3 = i2 + l;
                        im1[i3] /= sqrt_s1;
                        im1[i3 + 1] /= sqrt_s1;
                        im2[i3] /= sqrt_s2;
                        im2[i3 + 1] /= sqrt_s2;
                        l += 2;
                }
        }
        float * data = get_data();

        switch (mode) {
        case 0:
                for (int m1 = 0; m1 < 2; m1++) {
                        for (int m2 = 0; m2 < 2; m2++) {

                                if (m1 == 0 && m2 == 0) {
                                        for (int i = 0; i < steps; i++) {
                                                int i2 = i * rmax * 2;
                                                for (int j = 0; j < steps; j++) {
                                                        int l = i + j * steps * 2;
                                                        int j2 = j * rmax * 2;
                                                        data[l] = 0;
                                                        for (int k = 0; k < jmax; k++) {
                                                                data[l] += im1[i2 + k] * im2[j2 + k];
                                                        }
                                                }
                                        }
                                }
                                else {
                                        int steps2 = steps * m2 + steps * steps * 2 * m1;

                                        for (int i = 0; i < steps; i++) {
                                                int i2 = i * rmax * 2;
                                                for (int j = 0; j < steps; j++) {
                                                        int j2 = j * rmax * 2;
                                                        int l = i + j * steps * 2 + steps2;
                                                        data[l] = 0;

                                                        for (int k = 0; k < jmax; k += 2) {
                                                                i2 += k;
                                                                j2 += k;
                                                                data[l] += im1[i2] * im2[j2];
                                                                data[l] += -im1[i2 + 1] * im2[j2 + 1];
                                                        }
                                                }
                                        }
                                }
                        }
                }

                break;
        case 1:
                for (int m1 = 0; m1 < 2; m1++) {
                        for (int m2 = 0; m2 < 2; m2++) {
                                int steps2 = steps * m2 + steps * steps * 2 * m1;
                                int p1_sign = 1;
                                if (m1 != m2) {
                                        p1_sign = -1;
                                }

                                for (int i = 0; i < steps; i++) {
                                        int i2 = i * rmax * 2;

                                        for (int j = 0; j < steps; j++) {
                                                int j2 = j * rmax * 2;

                                                int l = i + j * steps * 2 + steps2;
                                                data[l] = 0;
                                                float a = 0;

                                                for (int k = 0; k < jmax; k += 2) {
                                                        i2 += k;
                                                        j2 += k;

#ifdef  _WIN32
                                                        float a1 = (float) _hypot(im1[i2], im1[i2 + 1]);
#else
                                                        float a1 = (float) hypot(im1[i2], im1[i2 + 1]);
#endif  //_WIN32
                                                        float p1 = atan2(im1[i2 + 1], im1[i2]);
                                                        float p2 = atan2(im2[j2 + 1], im2[j2]);

                                                        data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1;
                                                        a += a1;
                                                }

                                                data[l] /= (float)(a * M_PI / 180.0f);
                                        }
                                }
                        }
                }

                break;
        case 2:
                for (int m1 = 0; m1 < 2; m1++) {
                        for (int m2 = 0; m2 < 2; m2++) {
                                int steps2 = steps * m2 + steps * steps * 2 * m1;

                                for (int i = 0; i < steps; i++) {
                                        int i2 = i * rmax * 2;

                                        for (int j = 0; j < steps; j++) {
                                                int j2 = j * rmax * 2;
                                                int l = i + j * steps * 2 + steps2;
                                                data[l] = 0;

                                                for (int k = 0; k < jmax; k += 2) {
                                                        i2 += k;
                                                        j2 += k;
#ifdef  _WIN32
                                                        data[l] += (float) (_hypot(im1[i2], im1[i2 + 1]) * _hypot(im2[j2], im2[j2 + 1]));
#else
                                                        data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1]));
#endif  //_WIN32
                                                }
                                        }
                                }
                        }
                }

                break;
        default:
                break;
        }

        if (horizontal) {
                float *tmp_array = new float[ny];
                for (int i = 1; i < nx; i++) {
                        for (int j = 0; j < ny; j++) {
                                tmp_array[j] = get_value_at(i, j);
                        }
                        for (int j = 0; j < ny; j++) {
                                set_value_at(i, j, 0, tmp_array[(j + i) % ny]);
                        }
                }
                if( tmp_array )
                {
                        delete[]tmp_array;
                        tmp_array = 0;
                }
        }

        if( im1 )
        {
                delete[]im1;
                im1 = 0;
        }

        if( im2 )
        {
                delete im2;
                im2 = 0;
        }


        image1->update();
        image2->update();
        if( image1 )
        {
                delete image1;
                image1 = 0;
        }
        if( image2 )
        {
                delete image2;
                image2 = 0;
        }
        update();
        EXITFUNC;
}
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 3753 of file emdata.cpp.

References EMAN::EMData::copy(), data, ENTERFUNC, EXITFUNC, EMAN::EMData::get_data(), EMAN::EMData::get_ysize(), ImageFormatException, images, EMAN::EMUtil::is_same_size(), mean(), NullPointerException, EMAN::EMData::set_size(), EMAN::EMData::sqrt(), t, EMAN::EMData::transform(), and EMAN::EMData::update().

Referenced by main().

{
        ENTERFUNC;

        if (!image1 || !image2) {
                throw NullPointerException("NULL image");
        }

        if (!EMUtil::is_same_size(image1, image2)) {
                throw ImageFormatException("images not same size");
        }

        int steps2 = steps * 2;
        int image_ny = image1->get_ysize();
        EMData *image1_copy = image1->copy();
        EMData *image2_copy = image2->copy();

        float *im1 = new float[steps2 * image_ny];
        float *im2 = new float[steps2 * image_ny];

        EMData *images[] = { image1_copy, image2_copy };
        float *ims[] = { im1, im2 };

        for (int m = 0; m < 2; m++) {
                float *im = ims[m];
                float a = M_PI / steps2;
                Transform t(Dict("type","2d","alpha",-a));
                for (int i = 0; i < steps2; i++) {
                        images[i]->transform(t);
                        float *data = images[i]->get_data();

                        for (int j = 0; j < image_ny; j++) {
                                float sum = 0;
                                for (int k = 0; k < image_ny; k++) {
                                        sum += data[j * image_ny + k];
                                }
                                im[i * image_ny + j] = sum;
                        }

                        float sum1 = 0;
                        float sum2 = 0;
                        for (int j = 0; j < image_ny; j++) {
                                int l = i * image_ny + j;
                                sum1 += im[l];
                                sum2 += im[l] * im[l];
                        }

                        float mean = sum1 / image_ny;
                        float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1);

                        for (int j = 0; j < image_ny; j++) {
                                int l = i * image_ny + j;
                                im[l] = (im[l] - mean) / sigma;
                        }

                        images[i]->update();
                        a += M_PI / steps;
                }
        }

        set_size(steps2, steps2, 1);
        float *data1 = get_data();

        if (horiz) {
                for (int i = 0; i < steps2; i++) {
                        for (int j = 0, l = i; j < steps2; j++, l++) {
                                if (l == steps2) {
                                        l = 0;
                                }

                                float sum = 0;
                                for (int k = 0; k < image_ny; k++) {
                                        sum += im1[i * image_ny + k] * im2[l * image_ny + k];
                                }
                                data1[i + j * steps2] = sum;
                        }
                }
        }
        else {
                for (int i = 0; i < steps2; i++) {
                        for (int j = 0; j < steps2; j++) {
                                float sum = 0;
                                for (int k = 0; k < image_ny; k++) {
                                        sum += im1[i * image_ny + k] * im2[j * image_ny + k];
                                }
                                data1[i + j * steps2] = sum;
                        }
                }
        }

        update();

        if( image1_copy )
        {
                delete image1_copy;
                image1_copy = 0;
        }

        if( image2_copy )
        {
                delete image2_copy;
                image2_copy = 0;
        }

        if( im1 )
        {
                delete[]im1;
                im1 = 0;
        }

        if( im2 )
        {
                delete[]im2;
                im2 = 0;
        }
        EXITFUNC;
}