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