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. |
|
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.
Definition at line 3213 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(). 03215 { 03216 ENTERFUNC; 03217 03218 if (!image1 || !image2) { 03219 throw NullPointerException("NULL image"); 03220 } 03221 03222 if (mode < 0 || mode > 2) { 03223 throw OutofRangeException(0, 2, mode, "invalid mode"); 03224 } 03225 03226 if (!image1->is_complex()) { 03227 image1 = image1->do_fft(); 03228 } 03229 if (!image2->is_complex()) { 03230 image2 = image2->do_fft(); 03231 } 03232 03233 image1->ap2ri(); 03234 image2->ap2ri(); 03235 03236 if (!EMUtil::is_same_size(image1, image2)) { 03237 throw ImageFormatException("images not same sizes"); 03238 } 03239 03240 int image2_nx = image2->get_xsize(); 03241 int image2_ny = image2->get_ysize(); 03242 03243 int rmax = image2_ny / 4 - 1; 03244 int array_size = steps * rmax * 2; 03245 float *im1 = new float[array_size]; 03246 float *im2 = new float[array_size]; 03247 for (int i = 0; i < array_size; i++) { 03248 im1[i] = 0; 03249 im2[i] = 0; 03250 } 03251 03252 set_size(steps * 2, steps * 2, 1); 03253 03254 float *image1_data = image1->get_data(); 03255 float *image2_data = image2->get_data(); 03256 03257 float da = M_PI / steps; 03258 float a = -M_PI / 2.0f + da / 2.0f; 03259 int jmax = 0; 03260 03261 for (int i = 0; i < steps * 2; i += 2, a += da) { 03262 float s1 = 0; 03263 float s2 = 0; 03264 int i2 = i * rmax; 03265 int j = 0; 03266 03267 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) { 03268 float x = r * cos(a); 03269 float y = r * sin(a); 03270 03271 if (x < 0) { 03272 x = -x; 03273 y = -y; 03274 LOGERR("CCL ERROR %d, %f !\n", i, -x); 03275 } 03276 03277 int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx); 03278 int l = i2 + j; 03279 float x2 = x - floor(x); 03280 float y2 = y - floor(y); 03281 03282 im1[l] = Util::bilinear_interpolate(image1_data[k], 03283 image1_data[k + 2], 03284 image1_data[k + image2_nx], 03285 image1_data[k + 2 + image2_nx], x2, y2); 03286 03287 im2[l] = Util::bilinear_interpolate(image2_data[k], 03288 image2_data[k + 2], 03289 image2_data[k + image2_nx], 03290 image2_data[k + 2 + image2_nx], x2, y2); 03291 03292 k++; 03293 03294 im1[l + 1] = Util::bilinear_interpolate(image1_data[k], 03295 image1_data[k + 2], 03296 image1_data[k + image2_nx], 03297 image1_data[k + 2 + image2_nx], x2, y2); 03298 03299 im2[l + 1] = Util::bilinear_interpolate(image2_data[k], 03300 image2_data[k + 2], 03301 image2_data[k + image2_nx], 03302 image2_data[k + 2 + image2_nx], x2, y2); 03303 03304 s1 += Util::square_sum(im1[l], im1[l + 1]); 03305 s2 += Util::square_sum(im2[l], im2[l + 1]); 03306 } 03307 03308 jmax = j - 1; 03309 float sqrt_s1 = std::sqrt(s1); 03310 float sqrt_s2 = std::sqrt(s2); 03311 03312 int l = 0; 03313 for (float r = 1; r < rmax; r += 1.0f) { 03314 int i3 = i2 + l; 03315 im1[i3] /= sqrt_s1; 03316 im1[i3 + 1] /= sqrt_s1; 03317 im2[i3] /= sqrt_s2; 03318 im2[i3 + 1] /= sqrt_s2; 03319 l += 2; 03320 } 03321 } 03322 float * data = get_data(); 03323 03324 switch (mode) { 03325 case 0: 03326 for (int m1 = 0; m1 < 2; m1++) { 03327 for (int m2 = 0; m2 < 2; m2++) { 03328 03329 if (m1 == 0 && m2 == 0) { 03330 for (int i = 0; i < steps; i++) { 03331 int i2 = i * rmax * 2; 03332 for (int j = 0; j < steps; j++) { 03333 int l = i + j * steps * 2; 03334 int j2 = j * rmax * 2; 03335 data[l] = 0; 03336 for (int k = 0; k < jmax; k++) { 03337 data[l] += im1[i2 + k] * im2[j2 + k]; 03338 } 03339 } 03340 } 03341 } 03342 else { 03343 int steps2 = steps * m2 + steps * steps * 2 * m1; 03344 03345 for (int i = 0; i < steps; i++) { 03346 int i2 = i * rmax * 2; 03347 for (int j = 0; j < steps; j++) { 03348 int j2 = j * rmax * 2; 03349 int l = i + j * steps * 2 + steps2; 03350 data[l] = 0; 03351 03352 for (int k = 0; k < jmax; k += 2) { 03353 i2 += k; 03354 j2 += k; 03355 data[l] += im1[i2] * im2[j2]; 03356 data[l] += -im1[i2 + 1] * im2[j2 + 1]; 03357 } 03358 } 03359 } 03360 } 03361 } 03362 } 03363 03364 break; 03365 case 1: 03366 for (int m1 = 0; m1 < 2; m1++) { 03367 for (int m2 = 0; m2 < 2; m2++) { 03368 int steps2 = steps * m2 + steps * steps * 2 * m1; 03369 int p1_sign = 1; 03370 if (m1 != m2) { 03371 p1_sign = -1; 03372 } 03373 03374 for (int i = 0; i < steps; i++) { 03375 int i2 = i * rmax * 2; 03376 03377 for (int j = 0; j < steps; j++) { 03378 int j2 = j * rmax * 2; 03379 03380 int l = i + j * steps * 2 + steps2; 03381 data[l] = 0; 03382 float a = 0; 03383 03384 for (int k = 0; k < jmax; k += 2) { 03385 i2 += k; 03386 j2 += k; 03387 03388 #ifdef _WIN32 03389 float a1 = (float) _hypot(im1[i2], im1[i2 + 1]); 03390 #else 03391 float a1 = (float) hypot(im1[i2], im1[i2 + 1]); 03392 #endif //_WIN32 03393 float p1 = atan2(im1[i2 + 1], im1[i2]); 03394 float p2 = atan2(im2[j2 + 1], im2[j2]); 03395 03396 data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1; 03397 a += a1; 03398 } 03399 03400 data[l] /= (float)(a * M_PI / 180.0f); 03401 } 03402 } 03403 } 03404 } 03405 03406 break; 03407 case 2: 03408 for (int m1 = 0; m1 < 2; m1++) { 03409 for (int m2 = 0; m2 < 2; m2++) { 03410 int steps2 = steps * m2 + steps * steps * 2 * m1; 03411 03412 for (int i = 0; i < steps; i++) { 03413 int i2 = i * rmax * 2; 03414 03415 for (int j = 0; j < steps; j++) { 03416 int j2 = j * rmax * 2; 03417 int l = i + j * steps * 2 + steps2; 03418 data[l] = 0; 03419 03420 for (int k = 0; k < jmax; k += 2) { 03421 i2 += k; 03422 j2 += k; 03423 #ifdef _WIN32 03424 data[l] += (float) (_hypot(im1[i2], im1[i2 + 1]) * _hypot(im2[j2], im2[j2 + 1])); 03425 #else 03426 data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1])); 03427 #endif //_WIN32 03428 } 03429 } 03430 } 03431 } 03432 } 03433 03434 break; 03435 default: 03436 break; 03437 } 03438 03439 if (horizontal) { 03440 float *tmp_array = new float[ny]; 03441 for (int i = 1; i < nx; i++) { 03442 for (int j = 0; j < ny; j++) { 03443 tmp_array[j] = get_value_at(i, j); 03444 } 03445 for (int j = 0; j < ny; j++) { 03446 set_value_at(i, j, 0, tmp_array[(j + i) % ny]); 03447 } 03448 } 03449 if( tmp_array ) 03450 { 03451 delete[]tmp_array; 03452 tmp_array = 0; 03453 } 03454 } 03455 03456 if( im1 ) 03457 { 03458 delete[]im1; 03459 im1 = 0; 03460 } 03461 03462 if( im2 ) 03463 { 03464 delete im2; 03465 im2 = 0; 03466 } 03467 03468 03469 image1->update(); 03470 image2->update(); 03471 if( image1 ) 03472 { 03473 delete image1; 03474 image1 = 0; 03475 } 03476 if( image2 ) 03477 { 03478 delete image2; 03479 image2 = 0; 03480 } 03481 update(); 03482 EXITFUNC; 03483 }
|
|
Finds common lines between 2 real images.
Definition at line 3487 of file emdata.cpp. References EMAN::EMData::copy(), data, ENTERFUNC, EXITFUNC, EMAN::EMData::get_data(), ImageFormatException, images, EMAN::EMUtil::is_same_size(), NullPointerException, EMAN::EMData::set_size(), sqrt(), t, EMAN::EMData::transform(), and EMAN::EMData::update(). Referenced by main(). 03489 { 03490 ENTERFUNC; 03491 03492 if (!image1 || !image2) { 03493 throw NullPointerException("NULL image"); 03494 } 03495 03496 if (!EMUtil::is_same_size(image1, image2)) { 03497 throw ImageFormatException("images not same size"); 03498 } 03499 03500 int steps2 = steps * 2; 03501 int image_ny = image1->get_ysize(); 03502 EMData *image1_copy = image1->copy(); 03503 EMData *image2_copy = image2->copy(); 03504 03505 float *im1 = new float[steps2 * image_ny]; 03506 float *im2 = new float[steps2 * image_ny]; 03507 03508 EMData *images[] = { image1_copy, image2_copy }; 03509 float *ims[] = { im1, im2 }; 03510 03511 for (int m = 0; m < 2; m++) { 03512 float *im = ims[m]; 03513 float a = M_PI / steps2; 03514 Transform t(Dict("type","2d","alpha",-a)); 03515 for (int i = 0; i < steps2; i++) { 03516 images[i]->transform(t); 03517 float *data = images[i]->get_data(); 03518 03519 for (int j = 0; j < image_ny; j++) { 03520 float sum = 0; 03521 for (int k = 0; k < image_ny; k++) { 03522 sum += data[j * image_ny + k]; 03523 } 03524 im[i * image_ny + j] = sum; 03525 } 03526 03527 float sum1 = 0; 03528 float sum2 = 0; 03529 for (int j = 0; j < image_ny; j++) { 03530 int l = i * image_ny + j; 03531 sum1 += im[l]; 03532 sum2 += im[l] * im[l]; 03533 } 03534 03535 float mean = sum1 / image_ny; 03536 float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1); 03537 03538 for (int j = 0; j < image_ny; j++) { 03539 int l = i * image_ny + j; 03540 im[l] = (im[l] - mean) / sigma; 03541 } 03542 03543 images[i]->update(); 03544 a += M_PI / steps; 03545 } 03546 } 03547 03548 set_size(steps2, steps2, 1); 03549 float *data1 = get_data(); 03550 03551 if (horiz) { 03552 for (int i = 0; i < steps2; i++) { 03553 for (int j = 0, l = i; j < steps2; j++, l++) { 03554 if (l == steps2) { 03555 l = 0; 03556 } 03557 03558 float sum = 0; 03559 for (int k = 0; k < image_ny; k++) { 03560 sum += im1[i * image_ny + k] * im2[l * image_ny + k]; 03561 } 03562 data1[i + j * steps2] = sum; 03563 } 03564 } 03565 } 03566 else { 03567 for (int i = 0; i < steps2; i++) { 03568 for (int j = 0; j < steps2; j++) { 03569 float sum = 0; 03570 for (int k = 0; k < image_ny; k++) { 03571 sum += im1[i * image_ny + k] * im2[j * image_ny + k]; 03572 } 03573 data1[i + j * steps2] = sum; 03574 } 03575 } 03576 } 03577 03578 update(); 03579 03580 if( image1_copy ) 03581 { 03582 delete image1_copy; 03583 image1_copy = 0; 03584 } 03585 03586 if( image2_copy ) 03587 { 03588 delete image2_copy; 03589 image2_copy = 0; 03590 } 03591 03592 if( im1 ) 03593 { 03594 delete[]im1; 03595 im1 = 0; 03596 } 03597 03598 if( im2 ) 03599 { 03600 delete[]im2; 03601 im2 = 0; 03602 } 03603 EXITFUNC; 03604 }
|