casa5
Public
Commits
R. V. Urvashi authored c590664c200 Merge
| 310 310 | ++count; |
| 311 311 | } |
| 312 312 | } |
| 313 313 | } |
| 314 314 | } |
| 315 315 | } |
| 316 316 | |
| 317 317 | template<class T> SPIIT ImageRegridder<T>::_regridByVelocity() const { |
| 318 318 | const auto csysTo = this->_getTemplateCoords(); |
| 319 319 | const auto specCoordTo = csysTo.spectralCoordinate(); |
| 320 - | const auto specCoordFrom = this->_getImage()->coordinates() |
| 321 - | .spectralCoordinate(); |
| 320 + | const auto specCoordFrom |
| 321 + | = this->_getImage()->coordinates().spectralCoordinate(); |
| 322 322 | ThrowIf( |
| 323 323 | specCoordTo.frequencySystem(true) |
| 324 324 | != specCoordFrom.frequencySystem(true), |
| 325 325 | "Image to be regridded has different frequency system from template " |
| 326 326 | "coordinate system." |
| 327 327 | ); |
| 328 328 | ThrowIf( |
| 329 329 | specCoordTo.restFrequency() == 0, |
| 330 330 | "Template spectral coordinate rest frequency is 0, " |
| 331 331 | "so cannot regrid by velocity." |
| 332 332 | ); |
| 333 333 | ThrowIf( |
| 334 334 | specCoordFrom.restFrequency() == 0, |
| 335 335 | "Input image spectral coordinate rest frequency is 0, so cannot regrid " |
| 336 336 | "by velocity." |
| 337 337 | ); |
| 338 - | |
| 339 338 | std::unique_ptr<casacore::CoordinateSystem> csys( |
| 340 339 | dynamic_cast<casacore::CoordinateSystem *>(csysTo.clone()) |
| 341 340 | ); |
| 342 341 | auto templateSpecCoord = csys->spectralCoordinate(); |
| 343 342 | std::unique_ptr<casacore::CoordinateSystem> coordClone( |
| 344 343 | dynamic_cast<casacore::CoordinateSystem *>( |
| 345 344 | _subimage->coordinates().clone() |
| 346 345 | ) |
| 347 346 | ); |
| 348 347 | auto newSpecCoord = coordClone->spectralCoordinate(); |
| 349 348 | casacore::Double newVelRefVal = 0; |
| 350 - | casacore::Double newVelInc = 0; |
| 351 349 | std::pair<casacore::Double, casacore::Double> toVelLimits; |
| 352 350 | auto inSpecAxis = coordClone->spectralAxisNumber(false); |
| 353 - | for (uInt i=0; i<2; ++i) { |
| 351 + | casacore::Double newVelInc = 0.0; |
| 352 + | for (casacore::uInt i=0; i<2; ++i) { |
| 354 353 | // i == 0 => csysTo, i == 1 => csysFrom |
| 355 354 | auto *cs = i == 0 ? csys.get() : coordClone.get(); |
| 356 355 | // create and replace the coordinate system's spectral coordinate with |
| 357 356 | // a linear coordinate which describes the velocity axis. In this way |
| 358 357 | // we can regrid by velocity. |
| 359 358 | Int specCoordNum = cs->spectralCoordinateNumber(); |
| 360 359 | auto specCoord = cs->spectralCoordinate(); |
| 361 360 | if ( |
| 362 361 | specCoord.frequencySystem(false) != specCoord.frequencySystem(true) |
| 363 362 | ) { |
| 381 380 | if (cs == coordClone.get()) { |
| 382 381 | newSpecCoord = specCoord; |
| 383 382 | } |
| 384 383 | } |
| 385 384 | casacore::Double freqRefVal = specCoord.referenceValue()[0]; |
| 386 385 | casacore::Double velRefVal; |
| 387 386 | ThrowIf( |
| 388 387 | ! specCoord.frequencyToVelocity(velRefVal, freqRefVal), |
| 389 388 | "Unable to determine reference velocity" |
| 390 389 | ); |
| 391 - | casacore::Double vel0; |
| 390 + | casacore::Double vel0 = 0; |
| 392 391 | casacore::Double vel1 = 0; |
| 393 392 | ThrowIf( |
| 394 393 | ! specCoord.pixelToVelocity(vel0, 0.0) |
| 395 394 | || ! specCoord.pixelToVelocity(vel1, 1.0), |
| 396 395 | "Unable to determine velocity increment" |
| 397 396 | ); |
| 398 397 | if (i == 0) { |
| 399 398 | toVelLimits.first = vel0; |
| 400 399 | specCoord.pixelToVelocity( |
| 401 400 | toVelLimits.second, this->_getShape()[inSpecAxis] - 1 |
| 407 406 | if (i == 1) { |
| 408 407 | std::pair<casacore::Double, casacore::Double> fromVelLimits; |
| 409 408 | specCoord.pixelToVelocity(fromVelLimits.first, 0); |
| 410 409 | specCoord.pixelToVelocity( |
| 411 410 | fromVelLimits.second, _subimage->shape()[inSpecAxis] - 1 |
| 412 411 | ); |
| 413 412 | if (fromVelLimits.first > fromVelLimits.second) { |
| 414 413 | std::swap(fromVelLimits.first, fromVelLimits.second); |
| 415 414 | } |
| 416 415 | ThrowIf( |
| 417 - | fromVelLimits.first > toVelLimits.second |
| 418 - | || fromVelLimits.second < toVelLimits.first, |
| 416 + | ( |
| 417 + | fromVelLimits.first > toVelLimits.second |
| 418 + | && ! casacore::near(fromVelLimits.first, toVelLimits.second) |
| 419 + | ) |
| 420 + | || ( |
| 421 + | fromVelLimits.second < toVelLimits.first |
| 422 + | && ! casacore::near(fromVelLimits.second, toVelLimits.first) |
| 423 + | ), |
| 419 424 | "Request to regrid by velocity, but input and output velocity " |
| 420 425 | "coordinates do not overlap" |
| 421 426 | ); |
| 422 427 | } |
| 423 428 | casacore::Matrix<casacore::Double> pc(1, 1, 0); |
| 424 429 | pc.diagonal() = 1.0; |
| 425 430 | casacore::LinearCoordinate lin( |
| 426 431 | casacore::Vector<casacore::String>(1, "velocity"), |
| 427 432 | specCoord.worldAxisUnits(), |
| 428 433 | casacore::Vector<casacore::Double>(1, velRefVal), |
| 569 574 | auto nChan1 = md1.nChannels(); |
| 570 575 | casacore::Double world; |
| 571 576 | sp0.toWorld(world, 0); |
| 572 577 | auto end00 = world; |
| 573 578 | sp0.toWorld(world, nChan0 - 1); |
| 574 579 | auto end01 = world; |
| 575 580 | sp1.toWorld(world, 0); |
| 576 581 | auto end10 = world; |
| 577 582 | sp1.toWorld(world, nChan1 - 1); |
| 578 583 | auto end11 = world; |
| 584 + | auto minmax0 = minmax(end00, end01); |
| 585 + | auto minmax1 = minmax(end10, end11); |
| 579 586 | if ( |
| 580 - | max(end00, end01) < min(end10, end11) |
| 581 - | || max(end10, end11) < min(end00, end01) |
| 587 + | ( |
| 588 + | minmax0.second < minmax1.first |
| 589 + | && ! casacore::near(minmax0.second, minmax1.first) |
| 590 + | ) |
| 591 + | || ( |
| 592 + | minmax1.second < minmax0.first |
| 593 + | && ! casacore::near(minmax1.second, minmax0.first) |
| 594 + | ) |
| 582 595 | ) { |
| 583 596 | return false; |
| 584 597 | } |
| 585 598 | } |
| 586 599 | return true; |
| 587 600 | } |
| 588 601 | |
| 589 602 | template<class T> |
| 590 603 | casacore::Vector<std::pair<casacore::Double, casacore::Double>> |
| 591 604 | ImageRegridder<T>::_getDirectionCorners( |