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
	)			[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.

Parameters:

	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.

Exceptions:

	NullPointerException	If 'image1' or 'image2' is NULL.
	OutofRangeException	If 'mode' is invalid.
	ImageFormatException	If 'image1' 'image2' are not same size.

Definition at line 3197 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().

{
        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
	)			[inherited]

Finds common lines between 2 real images.

Parameters:

	image1	The first image.
	image2	The second image.
	steps	1/2 of the resolution of the map.
	horizontal	In horizontal way or not.

Exceptions:

	NullPointerException	If 'image1' or 'image2' is NULL.
	ImageFormatException	If 'image1' 'image2' are not same size.

Definition at line 3471 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(), NullPointerException, EMAN::EMData::set_size(), 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;
}