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

03199 {
03200         ENTERFUNC;
03201 
03202         if (!image1 || !image2) {
03203                 throw NullPointerException("NULL image");
03204         }
03205 
03206         if (mode < 0 || mode > 2) {
03207                 throw OutofRangeException(0, 2, mode, "invalid mode");
03208         }
03209 
03210         if (!image1->is_complex()) {
03211                 image1 = image1->do_fft();
03212         }
03213         if (!image2->is_complex()) {
03214                 image2 = image2->do_fft();
03215         }
03216 
03217         image1->ap2ri();
03218         image2->ap2ri();
03219 
03220         if (!EMUtil::is_same_size(image1, image2)) {
03221                 throw ImageFormatException("images not same sizes");
03222         }
03223 
03224         int image2_nx = image2->get_xsize();
03225         int image2_ny = image2->get_ysize();
03226 
03227         int rmax = image2_ny / 4 - 1;
03228         int array_size = steps * rmax * 2;
03229         float *im1 = new float[array_size];
03230         float *im2 = new float[array_size];
03231         for (int i = 0; i < array_size; i++) {
03232                 im1[i] = 0;
03233                 im2[i] = 0;
03234         }
03235 
03236         set_size(steps * 2, steps * 2, 1);
03237 
03238         float *image1_data = image1->get_data();
03239         float *image2_data = image2->get_data();
03240 
03241         float da = M_PI / steps;
03242         float a = -M_PI / 2.0f + da / 2.0f;
03243         int jmax = 0;
03244 
03245         for (int i = 0; i < steps * 2; i += 2, a += da) {
03246                 float s1 = 0;
03247                 float s2 = 0;
03248                 int i2 = i * rmax;
03249                 int j = 0;
03250 
03251                 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) {
03252                         float x = r * cos(a);
03253                         float y = r * sin(a);
03254 
03255                         if (x < 0) {
03256                                 x = -x;
03257                                 y = -y;
03258                                 LOGERR("CCL ERROR %d, %f !\n", i, -x);
03259                         }
03260 
03261                         int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx);
03262                         int l = i2 + j;
03263                         float x2 = x - floor(x);
03264                         float y2 = y - floor(y);
03265 
03266                         im1[l] = Util::bilinear_interpolate(image1_data[k],
03267                                                                                                 image1_data[k + 2],
03268                                                                                                 image1_data[k + image2_nx],
03269                                                                                                 image1_data[k + 2 + image2_nx], x2, y2);
03270 
03271                         im2[l] = Util::bilinear_interpolate(image2_data[k],
03272                                                                                                 image2_data[k + 2],
03273                                                                                                 image2_data[k + image2_nx],
03274                                                                                                 image2_data[k + 2 + image2_nx], x2, y2);
03275 
03276                         k++;
03277 
03278                         im1[l + 1] = Util::bilinear_interpolate(image1_data[k],
03279                                                                                                         image1_data[k + 2],
03280                                                                                                         image1_data[k + image2_nx],
03281                                                                                                         image1_data[k + 2 + image2_nx], x2, y2);
03282 
03283                         im2[l + 1] = Util::bilinear_interpolate(image2_data[k],
03284                                                                                                         image2_data[k + 2],
03285                                                                                                         image2_data[k + image2_nx],
03286                                                                                                         image2_data[k + 2 + image2_nx], x2, y2);
03287 
03288                         s1 += Util::square_sum(im1[l], im1[l + 1]);
03289                         s2 += Util::square_sum(im2[l], im2[l + 1]);
03290                 }
03291 
03292                 jmax = j - 1;
03293                 float sqrt_s1 = std::sqrt(s1);
03294                 float sqrt_s2 = std::sqrt(s2);
03295 
03296                 int l = 0;
03297                 for (float r = 1; r < rmax; r += 1.0f) {
03298                         int i3 = i2 + l;
03299                         im1[i3] /= sqrt_s1;
03300                         im1[i3 + 1] /= sqrt_s1;
03301                         im2[i3] /= sqrt_s2;
03302                         im2[i3 + 1] /= sqrt_s2;
03303                         l += 2;
03304                 }
03305         }
03306         float * data = get_data();
03307 
03308         switch (mode) {
03309         case 0:
03310                 for (int m1 = 0; m1 < 2; m1++) {
03311                         for (int m2 = 0; m2 < 2; m2++) {
03312 
03313                                 if (m1 == 0 && m2 == 0) {
03314                                         for (int i = 0; i < steps; i++) {
03315                                                 int i2 = i * rmax * 2;
03316                                                 for (int j = 0; j < steps; j++) {
03317                                                         int l = i + j * steps * 2;
03318                                                         int j2 = j * rmax * 2;
03319                                                         data[l] = 0;
03320                                                         for (int k = 0; k < jmax; k++) {
03321                                                                 data[l] += im1[i2 + k] * im2[j2 + k];
03322                                                         }
03323                                                 }
03324                                         }
03325                                 }
03326                                 else {
03327                                         int steps2 = steps * m2 + steps * steps * 2 * m1;
03328 
03329                                         for (int i = 0; i < steps; i++) {
03330                                                 int i2 = i * rmax * 2;
03331                                                 for (int j = 0; j < steps; j++) {
03332                                                         int j2 = j * rmax * 2;
03333                                                         int l = i + j * steps * 2 + steps2;
03334                                                         data[l] = 0;
03335 
03336                                                         for (int k = 0; k < jmax; k += 2) {
03337                                                                 i2 += k;
03338                                                                 j2 += k;
03339                                                                 data[l] += im1[i2] * im2[j2];
03340                                                                 data[l] += -im1[i2 + 1] * im2[j2 + 1];
03341                                                         }
03342                                                 }
03343                                         }
03344                                 }
03345                         }
03346                 }
03347 
03348                 break;
03349         case 1:
03350                 for (int m1 = 0; m1 < 2; m1++) {
03351                         for (int m2 = 0; m2 < 2; m2++) {
03352                                 int steps2 = steps * m2 + steps * steps * 2 * m1;
03353                                 int p1_sign = 1;
03354                                 if (m1 != m2) {
03355                                         p1_sign = -1;
03356                                 }
03357 
03358                                 for (int i = 0; i < steps; i++) {
03359                                         int i2 = i * rmax * 2;
03360 
03361                                         for (int j = 0; j < steps; j++) {
03362                                                 int j2 = j * rmax * 2;
03363 
03364                                                 int l = i + j * steps * 2 + steps2;
03365                                                 data[l] = 0;
03366                                                 float a = 0;
03367 
03368                                                 for (int k = 0; k < jmax; k += 2) {
03369                                                         i2 += k;
03370                                                         j2 += k;
03371 
03372 #ifdef  _WIN32
03373                                                         float a1 = (float) _hypot(im1[i2], im1[i2 + 1]);
03374 #else
03375                                                         float a1 = (float) hypot(im1[i2], im1[i2 + 1]);
03376 #endif  //_WIN32
03377                                                         float p1 = atan2(im1[i2 + 1], im1[i2]);
03378                                                         float p2 = atan2(im2[j2 + 1], im2[j2]);
03379 
03380                                                         data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1;
03381                                                         a += a1;
03382                                                 }
03383 
03384                                                 data[l] /= (float)(a * M_PI / 180.0f);
03385                                         }
03386                                 }
03387                         }
03388                 }
03389 
03390                 break;
03391         case 2:
03392                 for (int m1 = 0; m1 < 2; m1++) {
03393                         for (int m2 = 0; m2 < 2; m2++) {
03394                                 int steps2 = steps * m2 + steps * steps * 2 * m1;
03395 
03396                                 for (int i = 0; i < steps; i++) {
03397                                         int i2 = i * rmax * 2;
03398 
03399                                         for (int j = 0; j < steps; j++) {
03400                                                 int j2 = j * rmax * 2;
03401                                                 int l = i + j * steps * 2 + steps2;
03402                                                 data[l] = 0;
03403 
03404                                                 for (int k = 0; k < jmax; k += 2) {
03405                                                         i2 += k;
03406                                                         j2 += k;
03407 #ifdef  _WIN32
03408                                                         data[l] += (float) (_hypot(im1[i2], im1[i2 + 1]) * _hypot(im2[j2], im2[j2 + 1]));
03409 #else
03410                                                         data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1]));
03411 #endif  //_WIN32
03412                                                 }
03413                                         }
03414                                 }
03415                         }
03416                 }
03417 
03418                 break;
03419         default:
03420                 break;
03421         }
03422 
03423         if (horizontal) {
03424                 float *tmp_array = new float[ny];
03425                 for (int i = 1; i < nx; i++) {
03426                         for (int j = 0; j < ny; j++) {
03427                                 tmp_array[j] = get_value_at(i, j);
03428                         }
03429                         for (int j = 0; j < ny; j++) {
03430                                 set_value_at(i, j, 0, tmp_array[(j + i) % ny]);
03431                         }
03432                 }
03433                 if( tmp_array )
03434                 {
03435                         delete[]tmp_array;
03436                         tmp_array = 0;
03437                 }
03438         }
03439 
03440         if( im1 )
03441         {
03442                 delete[]im1;
03443                 im1 = 0;
03444         }
03445 
03446         if( im2 )
03447         {
03448                 delete im2;
03449                 im2 = 0;
03450         }
03451 
03452 
03453         image1->update();
03454         image2->update();
03455         if( image1 )
03456         {
03457                 delete image1;
03458                 image1 = 0;
03459         }
03460         if( image2 )
03461         {
03462                 delete image2;
03463                 image2 = 0;
03464         }
03465         update();
03466         EXITFUNC;
03467 }

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

03473 {
03474         ENTERFUNC;
03475 
03476         if (!image1 || !image2) {
03477                 throw NullPointerException("NULL image");
03478         }
03479 
03480         if (!EMUtil::is_same_size(image1, image2)) {
03481                 throw ImageFormatException("images not same size");
03482         }
03483 
03484         int steps2 = steps * 2;
03485         int image_ny = image1->get_ysize();
03486         EMData *image1_copy = image1->copy();
03487         EMData *image2_copy = image2->copy();
03488 
03489         float *im1 = new float[steps2 * image_ny];
03490         float *im2 = new float[steps2 * image_ny];
03491 
03492         EMData *images[] = { image1_copy, image2_copy };
03493         float *ims[] = { im1, im2 };
03494 
03495         for (int m = 0; m < 2; m++) {
03496                 float *im = ims[m];
03497                 float a = M_PI / steps2;
03498                 Transform t(Dict("type","2d","alpha",-a));
03499                 for (int i = 0; i < steps2; i++) {
03500                         images[i]->transform(t);
03501                         float *data = images[i]->get_data();
03502 
03503                         for (int j = 0; j < image_ny; j++) {
03504                                 float sum = 0;
03505                                 for (int k = 0; k < image_ny; k++) {
03506                                         sum += data[j * image_ny + k];
03507                                 }
03508                                 im[i * image_ny + j] = sum;
03509                         }
03510 
03511                         float sum1 = 0;
03512                         float sum2 = 0;
03513                         for (int j = 0; j < image_ny; j++) {
03514                                 int l = i * image_ny + j;
03515                                 sum1 += im[l];
03516                                 sum2 += im[l] * im[l];
03517                         }
03518 
03519                         float mean = sum1 / image_ny;
03520                         float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1);
03521 
03522                         for (int j = 0; j < image_ny; j++) {
03523                                 int l = i * image_ny + j;
03524                                 im[l] = (im[l] - mean) / sigma;
03525                         }
03526 
03527                         images[i]->update();
03528                         a += M_PI / steps;
03529                 }
03530         }
03531 
03532         set_size(steps2, steps2, 1);
03533         float *data1 = get_data();
03534 
03535         if (horiz) {
03536                 for (int i = 0; i < steps2; i++) {
03537                         for (int j = 0, l = i; j < steps2; j++, l++) {
03538                                 if (l == steps2) {
03539                                         l = 0;
03540                                 }
03541 
03542                                 float sum = 0;
03543                                 for (int k = 0; k < image_ny; k++) {
03544                                         sum += im1[i * image_ny + k] * im2[l * image_ny + k];
03545                                 }
03546                                 data1[i + j * steps2] = sum;
03547                         }
03548                 }
03549         }
03550         else {
03551                 for (int i = 0; i < steps2; i++) {
03552                         for (int j = 0; j < steps2; j++) {
03553                                 float sum = 0;
03554                                 for (int k = 0; k < image_ny; k++) {
03555                                         sum += im1[i * image_ny + k] * im2[j * image_ny + k];
03556                                 }
03557                                 data1[i + j * steps2] = sum;
03558                         }
03559                 }
03560         }
03561 
03562         update();
03563 
03564         if( image1_copy )
03565         {
03566                 delete image1_copy;
03567                 image1_copy = 0;
03568         }
03569 
03570         if( image2_copy )
03571         {
03572                 delete image2_copy;
03573                 image2_copy = 0;
03574         }
03575 
03576         if( im1 )
03577         {
03578                 delete[]im1;
03579                 im1 = 0;
03580         }
03581 
03582         if( im2 )
03583         {
03584                 delete[]im2;
03585                 im2 = 0;
03586         }
03587         EXITFUNC;
03588 }


Generated on Sun Mar 21 02:20:01 2010 for EMAN2 by  doxygen 1.5.6