31 #define M_PI 3.14159265358979323846 70 r =
sqrt(dx * dx + dy * dy);
75 return atan2(other->
dx * dy - dx * other->
dy, dx * other->
dx + dy * other->
dy);
79 stream <<
" dx " << segment.
dx <<
" dy " << segment.
dy <<
" r " << segment.
r <<
std::endl;
100 siStop = list.
begin();
106 for (
auto si1 = list.
begin(); si1 != siStop; ++si1, rank++)
107 for (
auto si2 = siStop; si2 != si1; --si2) {
108 push_back(
Segment(*si1, *si2, rank, transform));
110 sort(DecreasingLength);
130 for (SegmentPairListCIterator spi = pairList.
begin(); spi != pairList.
end(); spi++) {
132 if (a_pair.first->s1rank != rank1 || a_pair.second->s1rank != rank2)
continue;
136 if (matchList->size() == 0)
137 matchList->push_back(
StarMatch(transform.
apply(*(a_pair.first->s1)), *(a_pair.second->s1),
138 a_pair.first->s1, a_pair.second->s1));
140 matchList->push_back(
StarMatch(transform.
apply(*(a_pair.first->s2)), *(a_pair.second->s2),
141 a_pair.first->s2, a_pair.second->s2));
148 int idiff = first->size() - second->size();
152 return (first->getDist2() < second->getDist2());
174 int nBinsAngle = 180;
175 double angleOffset =
M_PI / nBinsAngle;
178 Histo2d histo(nBinsR, minRatio, maxRatio, nBinsAngle, -
M_PI - angleOffset,
M_PI - angleOffset);
180 SegmentIterator segi1, segi2;
183 for (segi1 = sList1.
begin(); segi1 != sList1.
end(); ++segi1) {
185 if (seg1->
r == 0)
continue;
186 for (segi2 = sList2.
begin(); segi2 != sList2.
end(); ++segi2) {
193 ratio = seg2->
r / seg1->
r;
194 if (ratio > maxRatio)
continue;
195 if (ratio < minRatio)
break;
197 if (angle >
M_PI - angleOffset) angle -= 2. *
M_PI;
198 histo.
fill(ratio, angle);
212 double ratioMax, angleMax;
213 double maxContent = histo.
maxBin(ratioMax, angleMax);
214 histo.
fill(ratioMax, angleMax, -maxContent);
217 LOGLS_DEBUG(_log,
" valMax " << maxContent <<
" ratio " << ratioMax <<
" angle " << angleMax);
219 minRatio = ratioMax - binr / 2;
220 maxRatio = ratioMax + binr / 2;
221 double minAngle = angleMax - bina / 2;
222 double maxAngle = angleMax + bina / 2;
223 SegmentPairList pairList;
228 for (segi1 = sList1.
begin(); segi1 != sList1.
end(); ++segi1) {
230 if (seg1->
r == 0)
continue;
231 for (segi2 = sList2.
begin(); segi2 != sList2.
end(); ++segi2) {
233 ratio = seg2->
r / seg1->
r;
234 if (ratio > maxRatio)
continue;
235 if (ratio < minRatio)
238 if (angle >
M_PI - angleOffset) angle -= 2. *
M_PI;
239 if (angle < minAngle || angle > maxAngle)
continue;
244 for (
int iteration = 0; iteration < conditions.
maxTrialCount; iteration++) {
246 double maxval = historank.
maxBin(dr1, dr2);
248 historank.
fill(dr1, dr2, -maxval);
250 a_list->refineTransform(conditions.
nSigmas);
254 Solutions.
sort(DecreasingQuality);
260 LOGLS_DEBUG(_log,
"Best solution " << best->computeResidual() <<
" npairs " << best->size());
262 LOGLS_DEBUG(_log,
"Chi2 " << best->getChi2() <<
',' <<
" Number of solutions " << Solutions.
size());
279 if (list1.
size() <= 4 || list2.
size() <= 4) {
280 LOGL_FATAL(_log,
"ListMatchupRotShift_New : (at least) one of the lists is too short.");
293 int nBinsAngle = 180;
294 double angleOffset =
M_PI / nBinsAngle;
301 SegmentIterator segi1, segi2;
305 for (segi1 = sList1.
begin(); segi1 != sList1.
end(); ++segi1) {
307 if (seg1->
r == 0)
continue;
308 for (segi2 = sList2.
begin(); segi2 != sList2.
end(); ++segi2) {
315 ratio = seg2->
r / seg1->
r;
316 if (ratio > maxRatio)
continue;
317 if (ratio < minRatio)
break;
319 if (angle >
M_PI - angleOffset) angle -= 2. *
M_PI;
320 histo.fill(ratio, angle, seg1->
s1rank + 0.5, seg2->
s1rank + 0.5);
331 int oldMaxContent = 0;
337 int maxContent = histo.maxBin(pars);
338 if (maxContent == 0)
break;
340 LOGLS_DEBUG(_log,
"ValMax " << maxContent <<
" ratio " << pars[0] <<
" angle " << pars[1]);
347 if (maxContent < oldMaxContent && i >= conditions.
maxTrialCount)
break;
349 oldMaxContent = maxContent;
351 int rank1L1 = int(pars[2]);
352 int rank1L2 = int(pars[3]);
353 double minAngle, maxAngle;
354 histo.binLimits(pars, 0, minRatio, maxRatio);
355 histo.binLimits(pars, 1, minAngle, maxAngle);
359 for (segi1 = sList1.
begin(); segi1 != sList1.
end(); ++segi1) {
361 if (seg1->
s1rank != rank1L1)
continue;
362 if (seg1->
r == 0)
continue;
363 for (segi2 = sList2.
begin(); segi2 != sList2.
end(); ++segi2) {
365 if (seg2->
s1rank != rank1L2)
continue;
367 if (a_list->size() == 0)
369 ratio = seg2->
r / seg1->
r;
370 if (ratio > maxRatio)
continue;
371 if (ratio < minRatio)
374 if (angle >
M_PI - angleOffset) angle -= 2. *
M_PI;
375 if (angle < minAngle || angle > maxAngle)
continue;
387 if (
int(a_list->size()) != maxContent + 1) {
388 LOGLS_ERROR(_log,
"There is an internal inconsistency in ListMatchupRotShift.");
390 LOGLS_ERROR(_log,
"matches->size() = " << a_list->size());
392 a_list->refineTransform(conditions.
nSigmas);
396 if (Solutions.
size() == 0) {
397 LOGLS_ERROR(_log,
"Error In ListMatchup : not a single pair match.");
398 LOGLS_ERROR(_log,
"Probably, the relative scale of lists is not within bounds.");
399 LOGLS_ERROR(_log,
"min/max ratios: " << minRatio <<
' ' << maxRatio);
403 Solutions.
sort(DecreasingQuality);
409 LOGLS_INFO(_log,
"Best solution " << best->computeResidual() <<
" npairs " << best->size());
411 LOGLS_INFO(_log,
"Chi2 " << best->getChi2() <<
", Number of solutions " << Solutions.
size());
420 return ListMatchupRotShift_Old(list1, list2, transform, conditions);
422 return ListMatchupRotShift_New(list1, list2, transform, conditions);
445 "unflipped Residual " << unflipped->computeResidual() <<
" nused " << unflipped->size());
446 LOGLS_DEBUG(_log,
"flipped Residual " << flipped->computeResidual() <<
" nused " << flipped->size());
448 if (DecreasingQuality(flipped, unflipped)) {
466 int ncomb = list1.
size() * list2.
size();
467 if (!ncomb)
return nullptr;
472 nx = (int)
sqrt(ncomb);
474 Histo2d histo(nx, -maxShift, maxShift, nx, -maxShift, maxShift);
478 for (s1 = list1.
begin(); s1 != list1.
end(); ++
s1) {
479 transform.
apply((*s1)->x, (*s1)->y, x1, y1);
480 for (s2 = list2.
begin(); s2 != list2.
end(); ++
s2) {
481 histo.
fill((*s2)->x - x1, (*s2)->y - y1);
484 double dx = 0,
dy = 0;
494 double maxShift,
double binSize) {
497 int ncomb = list1.
size() * list2.
size();
501 nx = (int)
sqrt(
double(ncomb));
504 nx = int(2 * maxShift / binSize + 0.5);
506 Histo2d histo(nx, -maxShift, maxShift, nx, -maxShift, maxShift);
507 double binSizeNew = 2 * maxShift / nx;
512 for (s1 = list1.
begin(); s1 != list1.
end(); ++
s1) {
513 transform.
apply((*s1)->x, (*s1)->y, x1, y1);
522 for (
int i = 0; i < 4; ++i) {
523 double dx = 0,
dy = 0;
525 histo.
fill(dx,
dy, -count);
527 auto newGuess =
compose(shift, transform);
528 auto raw_matches =
listMatchCollect(list1, list2, newGuess.get(), binSizeNew);
530 raw_matches->applyTransform(*matches, &transform);
531 matches->setTransformOrder(1);
532 matches->refineTransform(3.);
535 Solutions.
sort(DecreasingQuality);
538 std::dynamic_pointer_cast<const AstrometryTransformLinear>(
539 Solutions.
front()->getTransform()))));
552 const Point *p1 = (*si);
554 const BaseStar *neighbour = list2.findClosest(p2);
555 if (!neighbour)
continue;
557 if (distance < maxDist) {
558 matches->push_back(
StarMatch(*p1, *neighbour, *si, neighbour));
560 matches->back().distance =
distance;
578 if (!neighbour)
continue;
580 if (distance < maxDist) {
581 matches->push_back(
StarMatch(*p1, *neighbour, p1, neighbour));
583 matches->back().distance =
distance;
586 matches->setTransform(guess);
600 m->setTransform(bestTransform);
601 m->refineTransform(3.);
602 LOGLS_INFO(_log,
"Iterating: resid " <<
m->computeResidual() <<
" size " <<
m->size());
604 (prevMatch &&
m->computeResidual() < prevMatch->computeResidual() * 0.999 &&
m->Chi2() > 0)) {
606 bestTransform = prevMatch->Transform();
616 const double maxDist) {
622 if (!neighbour)
continue;
623 double distance = p1->Distance(*neighbour);
624 if (distance < maxDist) {
625 matches->push_back(
StarMatch(*p1, *neighbour, p1, neighbour));
627 matches->back().distance =
distance;
631 matches->setTransform(std::make_shared<AstrometryTransformIdentity>());
636 static bool is_transform_ok(
const StarMatchList *match,
double pixSizeRatio2,
const size_t nmin) {
637 if ((
fabs(
fabs(std::dynamic_pointer_cast<const AstrometryTransformLinear>(match->
getTransform())
642 (match->
size() > nmin))
660 double dx = tf1.
x - tf2.
x;
661 double dy = tf1.
y - tf2.
y;
662 diff2 += (tf1.
vy * dx * dx + tf1.
vx * dy * dy - 2 * tf1.
vxy * dx *
dy) /
666 if (count)
return diff2 / double(count);
671 size_t nstars = match->
size();
674 for (
auto it = match->
begin(); it != match->
end(); ++it, ++ir)
675 *ir =
sqrt(transform->
apply(it->point1).computeDist2(it->point2));
677 return (nstars & 1) ? resid[nstars / 2] : (resid[nstars / 2 - 1] + resid[nstars / 2]) * 0.5;
689 LOGLS_INFO(_log,
"listMatchCombinatorial: find match between " << list1.
size() <<
" and " << list2.
size()
697 if (is_transform_ok(match.get(), pixSizeRatio2, nmin))
698 transform = match->getTransform()->clone();
700 LOGL_ERROR(_log,
"listMatchCombinatorial: direct transform failed, trying reverse");
702 if (is_transform_ok(match.get(), pixSizeRatio2, nmin))
703 transform = match->inverseTransform();
712 LOGL_DEBUG(_log,
" listMatchCombinatorial: found the following transform.");
716 LOGL_ERROR(_log,
"listMatchCombinatorial: failed to find a transform");
722 const int maxOrder) {
728 const double brightDist = 2.;
729 const double fullDist = 4.;
730 const double nSigmas = 3.;
731 const size_t nStars = 500;
746 double curChi2 =
computeChi2(*brightMatch, *transform) / brightMatch->size();
748 LOGLS_INFO(_log,
"listMatchRefine: start: med.resid " << median_distance(fullMatch.get(), transform.
get())
749 <<
" #match " << fullMatch->size());
752 auto curTransform = brightMatch->getTransform()->clone();
756 brightMatch->setTransformOrder(order);
757 brightMatch->refineTransform(nSigmas);
758 transDiff = transform_diff(list1, brightMatch->getTransform().get(), curTransform.get());
759 curTransform = brightMatch->getTransform()->clone();
760 brightMatch =
listMatchCollect(list1, list2, curTransform.get(), brightDist);
761 }
while (brightMatch->size() > nstarmin && transDiff > 0.05 && ++iter < 5);
763 double prevChi2 = curChi2;
764 curChi2 =
computeChi2(*brightMatch, *curTransform) / brightMatch->size();
767 LOGLS_INFO(_log,
"listMatchRefine: order " << order <<
" med.resid " 768 << median_distance(fullMatch.get(), curTransform.get())
769 <<
" #match " << fullMatch->size());
770 if (((prevChi2 - curChi2) > 0.01 * curChi2) && curChi2 > 0) {
771 LOGLS_INFO(_log,
" listMatchRefine: order " << order <<
" was a better guess.");
772 transform = brightMatch->getTransform()->clone();
774 nstarmin = brightMatch->getTransform()->getNpar();
775 }
while (++order <= maxOrder);
BaseStarList::const_iterator BaseStarCIterator
std::shared_ptr< const BaseStar > s1
void binWidth(double &Hdx, double &Hdy) const
SegmentPair(Segment *f, Segment *s)
A hanger for star associations.
std::shared_ptr< const BaseStar > s2
std::list< Segment >::const_iterator SegmentCIterator
void fill(float x, float y, float weight=1.)
A class to histogram in 4 dimensions.
double maxSizeRatio() const
table::Key< double > angle
#define LOGL_ERROR(logger, message...)
Log a error-level message using a varargs/printf style interface.
ItemVariant const * other
std::unique_ptr< StarMatchList > listMatchCollect(const BaseStarList &list1, const BaseStarList &list2, const AstrometryTransform *guess, const double maxDist)
assembles star matches.
double Distance(const Point &other) const
SegmentPairList::const_iterator SegmentPairListCIterator
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
A Point with uncertainties.
void cutTail(const int nKeep)
cuts the end of the std::list
std::unique_ptr< AstrometryTransform > listMatchCombinatorial(const BaseStarList &list1, const BaseStarList &list2, const MatchConditions &conditions=MatchConditions())
std::unique_ptr< StarMatchList > matchSearchRotShift(BaseStarList &list1, BaseStarList &list2, const MatchConditions &conditions)
searches a geometrical transformation that goes from list1 to list2.
std::unique_ptr< AstrometryTransformLinear > listMatchupShift(const BaseStarList &list1, const BaseStarList &list2, const AstrometryTransform &transform, double maxShift, double binSize=0)
searches for a 2 dimensional shift using a very crude histogram method.
double maxBin(double &x, double &y) const
LSST DM logging module built on log4cxx.
std::unique_ptr< AstrometryTransform > listMatchRefine(const BaseStarList &list1, const BaseStarList &list2, std::unique_ptr< AstrometryTransform > transform, const int maxOrder=3)
The base class for handling stars. Used by all matching routines.
std::shared_ptr< const BaseStar > findClosest(const Point &where, const double maxDist, bool(*SkipIt)(const BaseStar &)=nullptr) const
Find the closest with some rejection capability.
std::unique_ptr< StarMatchList > matchSearchRotShiftFlip(BaseStarList &list1, BaseStarList &list2, const MatchConditions &conditions)
same as above but searches also a flipped solution.
#define LOGL_INFO(logger, message...)
Log a info-level message using a varargs/printf style interface.
#define LOGL_FATAL(logger, message...)
Log a fatal-level message using a varargs/printf style interface.
std::list< Segment >::iterator SegmentIterator
#define LOGLS_DEBUG(logger, message)
Log a debug-level message using an iostream-based interface.
A base class for image defects.
double relativeAngle(Segment *other)
Iterator beginScan(const Point &where, double maxDist) const
double minSizeRatio() const
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
std::list< SegmentPair > SegmentPairList
Iterator meant to traverse objects within some limiting distance.
void dumpTransform(std::ostream &stream=std::cout) const
print the matching transformation quality (transform, chi2, residual)
#define LOGLS_INFO(logger, message)
Log a info-level message using an iostream-based interface.
void copyTo(StarList< Star > ©) const
clears copy and makes a copy of the std::list to copy
std::shared_ptr< const AstrometryTransform > getTransform() const
carries out a fit with outlier rejection
void fluxSort()
a model routine to sort the std::list
This is an auxillary class for matching objects from starlists.
Combinatorial searches for linear transformations to go from list1 to list2.
Segment(std::shared_ptr< const BaseStar > star1, std::shared_ptr< const BaseStar > star2, const int star1Rank, const AstrometryTransform &transform)
friend std::ostream & operator<<(std::ostream &stream, const Segment &segment)
SegmentList(const BaseStarList &list, const int nStar, const AstrometryTransform &transform=AstrometryTransformIdentity())
#define LOG_GET(logger)
Returns a Log object associated with logger.
Parameters to be provided to combinatorial searches.
double computeChi2(const StarMatchList &L, const AstrometryTransform &transform)
the actual chi2
Fast locator in starlists.
#define LOGLS_ERROR(logger, message)
Log a error-level message using an iostream-based interface.
SegmentPairList::iterator SegmentPairListIterator