9 #include "boost/shared_array.hpp" 10 #include "boost/multi_index_container.hpp" 11 #include "boost/multi_index/sequenced_index.hpp" 12 #include "boost/multi_index/ordered_index.hpp" 13 #include "boost/multi_index/global_fun.hpp" 15 #include "gsl/gsl_linalg.h" 34 struct insertionOrderTag{};
35 typedef boost::multi_index_container<
37 boost::multi_index::indexed_by<
38 boost::multi_index::sequenced<
39 boost::multi_index::tag<insertionOrderTag>
41 boost::multi_index::ordered_unique<
42 boost::multi_index::tag<refIdTag>,
43 boost::multi_index::global_fun<
47 > MultiIndexedProxyPairList;
58 inline double absDeltaAngle(
double ang1,
double ang2) {
62 bool cmpPair(ProxyPair
const &
a, ProxyPair
const &
b) {
71 struct CompareProxyFlux {
74 double aFlux = a.
record->get(key);
75 double bFlux = b.
record->get(key);
104 if (startInd >= a.
size()) {
107 CompareProxyFlux cmp = {key};
111 return ProxyVector(b.begin() + startInd, b.begin() + endInd);
122 for (
size_t i = 0; i < a.
size(); i++) {
124 double dpa = absDeltaAngle(a[i].pa, p.
pa);
125 if (dd < e && dpa < e_dpa) {
142 double dd_min = 1.E+10;
145 for (
auto i = a.
begin(); i != a.
end(); ++i) {
149 absDeltaAngle(p.
pa, i->pa - dpa) < e_dpa &&
151 (i->first == q.
first)) {
157 absDeltaAngle(p.
pa, i->pa - dpa) < dpa_min) {
169 boost::shared_array<double>
const &
coeff,
175 int ncoeff = (order + 1) * (order + 2) / 2;
179 for (
int i = 0; i <= order; i++) {
180 for (
int k = 0; k <= i; k++) {
182 *xn += coeff[n] *
pow(x, j) *
pow(y, k);
189 boost::shared_array<double> polyfit(
194 int ncoeff = (order + 1) * (order + 2) / 2;
199 for (
int i = 0; i <= order; i++) {
200 for (
int k = 0; k <= i; k++) {
209 for (
size_t k = 0; k < img.
size(); k++) {
217 boost::shared_array<double>
coeff(
new double[ncoeff*2]);
219 for (
int loop = 0; loop < 1; loop++) {
220 for (
int i = 0; i <
ncoeff; i++) {
221 for (
int j = 0; j <
ncoeff; j++) {
222 a_data[i*ncoeff+j] = 0.0;
223 for (
size_t k = 0; k < img.
size(); k++) {
225 a_data[i*ncoeff+j] +=
pow(img[k].getX(), xorder[i]) *
226 pow(img[k].getY(), yorder[i]) *
227 pow(img[k].getX(), xorder[j]) *
228 pow(img[k].getY(), yorder[j]);
232 b_data[i] = c_data[i] = 0.0;
233 for (
unsigned int k = 0; k < img.
size(); k++) {
235 b_data[i] +=
pow(img[k].getX(), xorder[i]) *
236 pow(img[k].getY(), yorder[i]) *
238 c_data[i] +=
pow(img[k].getX(), xorder[i]) *
239 pow(img[k].getY(), yorder[i]) *
245 gsl_matrix_view a = gsl_matrix_view_array(a_data.get(),
ncoeff,
ncoeff);
246 gsl_vector_view b = gsl_vector_view_array(b_data.get(),
ncoeff);
247 gsl_vector_view c = gsl_vector_view_array(c_data.get(),
ncoeff);
256 gsl_linalg_LU_decomp(&a.matrix, p.get(), &
s);
257 gsl_linalg_LU_solve(&a.matrix, p.get(), &b.vector, x.get());
258 gsl_linalg_LU_solve(&a.matrix, p.get(), &c.vector, y.get());
260 for (
int i = 0; i <
ncoeff; i++) {
261 coeff[i] = x->data[i];
262 coeff[i+
ncoeff] = y->data[i];
265 double S, Sx, Sy, Sxx, Syy;
266 S = Sx = Sy = Sxx = Syy = 0.0;
267 for (
size_t k = 0; k < img.
size(); k++) {
269 double x0 = img[k].getX();
270 double y0 = img[k].getY();
272 transform(order, coeff, x0, y0, &x1, &y1);
274 Sx += (x1 - posRefCat[k].getX());
275 Sxx += (x1 - posRefCat[k].getX()) * (x1 - posRefCat[k].getX());
276 Sy += (y1 - posRefCat[k].getY());
277 Syy += (y1 - posRefCat[k].getY()) * (y1 - posRefCat[k].getY());
280 double x_sig =
std::sqrt((Sxx - Sx * Sx / S) / S);
281 double y_sig =
std::sqrt((Syy - Sy * Sy / S) / S);
284 for (
size_t k = 0; k < img.
size(); k++) {
285 double x0 = img[k].getX();
286 double y0 = img[k].getY();
288 transform(order, coeff, x0, y0, &x1, &y1);
289 if (
std::fabs(x1-posRefCat[k].getX()) > 2. * x_sig ||
290 std::fabs(y1-posRefCat[k].getY()) > 2. * y_sig) {
316 double matchingAllowancePix
318 auto minDistSq = matchingAllowancePix*matchingAllowancePix;
319 auto foundPtr = posRefCat.
end();
320 for (
auto posRefPtr = posRefCat.
begin(); posRefPtr != posRefCat.
end(); ++posRefPtr) {
321 auto const dx = posRefPtr->getX() -
x;
322 auto const dy = posRefPtr->getY() - y;
323 auto const distSq = dx*dx + dy*dy;
324 if (distSq < minDistSq) {
325 foundPtr = posRefPtr;
350 void addNearestMatch(
351 MultiIndexedProxyPairList &proxyPairList,
352 boost::shared_array<double> coeff,
355 double matchingAllowancePix
358 auto x0 = source.
getX();
359 auto y0 = source.
getY();
361 auto refObjDist = searchNearestPoint(posRefCat, x1, y1, matchingAllowancePix);
362 if (refObjDist.first == posRefCat.
end()) {
366 auto existingMatch = proxyPairList.get<refIdTag>().
find(refObjDist.first->record->getId());
367 if (existingMatch == proxyPairList.get<refIdTag>().end()) {
369 auto proxyPair = ProxyPair(*refObjDist.first, source);
370 proxyPairList.get<refIdTag>().insert(proxyPair);
378 if (existingMatch->distance <= refObjDist.second) {
383 proxyPairList.get<refIdTag>().
erase(existingMatch);
384 auto proxyPair = ProxyPair(*refObjDist.first, source);
385 proxyPairList.get<refIdTag>().insert(proxyPair);
401 boost::shared_array<double> coeff,
404 double matchingAllowancePix,
407 MultiIndexedProxyPairList proxyPairList;
409 for (
auto sourcePtr = sourceCat.
begin(); sourcePtr != sourceCat.
end(); ++sourcePtr) {
410 addNearestMatch(proxyPairList, coeff, posRefCat, *sourcePtr, matchingAllowancePix);
413 if (proxyPairList.size() > 5) {
414 for (
auto j = 0; j < 100; j++) {
415 auto prevNumMatches = proxyPairList.size();
418 srcMat.reserve(proxyPairList.size());
419 catMat.reserve(proxyPairList.size());
420 for (
auto matchPtr = proxyPairList.get<refIdTag>().begin();
421 matchPtr != proxyPairList.get<refIdTag>().
end(); ++matchPtr) {
422 catMat.push_back(matchPtr->first);
423 srcMat.push_back(matchPtr->second);
425 coeff = polyfit(order, srcMat, catMat);
426 proxyPairList.clear();
428 for (ProxyVector::const_iterator sourcePtr = sourceCat.
begin();
429 sourcePtr != sourceCat.
end(); ++sourcePtr) {
430 addNearestMatch(proxyPairList, coeff, posRefCat, *sourcePtr, matchingAllowancePix);
432 if (proxyPairList.size() == prevNumMatches) {
438 matPair.reserve(proxyPairList.size());
439 for (
auto proxyPairIter = proxyPairList.get<insertionOrderTag>().begin();
440 proxyPairIter != proxyPairList.get<insertionOrderTag>().
end(); ++proxyPairIter) {
442 proxyPairIter->first.record,
443 std::static_pointer_cast<afwTable::SourceRecord>(proxyPairIter->second.record),
444 proxyPairIter->distance
458 if (refFluxField.empty()) {
461 if (sourceFluxField.empty()) {
464 if (numBrightStars <= 0) {
467 if (minMatchedPairs < 0) {
470 if (matchingAllowancePix <= 0) {
473 if (maxOffsetPix <= 0) {
476 if (maxRotationDeg <= 0) {
479 if (allowedNonperpDeg <= 0) {
482 if (numPointsForShape <= 0) {
485 if (maxDeterminant <= 0) {
497 for (afwTable::SourceCatalog::const_iterator sourcePtr = sourceCat.begin();
498 sourcePtr != sourceCat.end(); ++sourcePtr) {
513 posRefPtr != posRefCat.
end(); ++posRefPtr) {
528 if (posRefCat.
empty()) {
531 if (sourceCat.empty()) {
534 if (posRefBegInd < 0) {
537 if (static_cast<size_t>(posRefBegInd) >= posRefCat.
size()) {
545 for (
auto iter = sourceCat.begin(); iter != sourceCat.end(); ++iter) {
548 srcCenter /= sourceCat.size();
551 for (
auto iter = posRefCat.
begin(); iter != posRefCat.
end(); ++iter) {
554 refCenter /= posRefCat.
size();
576 sourceSubCat.
size()+25,
584 size_t const posRefCatSubSize = posRefSubCat.
size();
585 for (
size_t i = 0; i < posRefCatSubSize-1; i++) {
586 for (
size_t j = i+1; j < posRefCatSubSize; j++) {
587 posRefPairList.
push_back(ProxyPair(posRefSubCat[i], posRefSubCat[j]));
594 size_t const sourceSubCatSize = sourceSubCat.
size();
595 for (
size_t i = 0; i < sourceSubCatSize-1; i++) {
596 for (
size_t j = i+1; j < sourceSubCatSize; j++) {
597 sourcePairList.push_back(ProxyPair(sourceSubCat[i], sourceSubCat[j]));
600 std::sort(sourcePairList.begin(), sourcePairList.end(), cmpPair);
607 for (
size_t ii = 0; ii < sourcePairList.size(); ii++) {
608 ProxyPair p = sourcePairList[ii];
620 for (
size_t l = 0; l < q.
size(); l++) {
623 double dpa = p.
pa - q[l].pa;
636 for (
size_t k = 0; k < sourceSubCat.
size(); k++) {
637 if (p.
first == sourceSubCat[k] || p.
second == sourceSubCat[k])
continue;
639 ProxyPair pp(p.
first, sourceSubCat[k]);
643 if (r != posRefPairList.
end()) {
650 if (srcMatPair.
size() == fullShapeSize) {
656 bool goodMatch =
false;
657 if (srcMatPair.
size() == fullShapeSize) {
659 for (
size_t k = 1; k < catMatPair.
size(); k++) {
660 if (catMatPair[0].
first != catMatPair[k].
first) {
667 if (goodMatch && srcMatPair.
size() == fullShapeSize) {
674 for (
size_t k = 0; k < srcMatPair.
size(); k++) {
679 boost::shared_array<double> coeff = polyfit(1, srcMat, catMat);
682 for (
size_t k = 0; k < srcMat.
size(); k++) {
683 std::cout <<
"circle(" << srcMat[k].getX() <<
"," 684 << srcMat[k].getY() <<
",10) # color=green" <<
std::endl;
685 std::cout <<
"circle(" << catMat[k].getX() <<
"," 686 << catMat[k].getY() <<
",10) # color=red" <<
std::endl;
687 std::cout <<
"line(" << srcMat[0].getX() <<
"," << srcMat[0].getY() <<
"," 688 << srcMat[k].getX() <<
"," << srcMat[k].getY()
689 <<
") # line=0 0 color=green" <<
std::endl;
690 std::cout <<
"line(" << catMat[0].getX() <<
"," << catMat[0].getY() <<
"," 691 << catMat[k].getX() <<
"," << catMat[k].getY()
692 <<
") # line=0 0 color=red" <<
std::endl;
709 std::cout <<
"Angle between axes (deg; allowed 90 +/- ";
713 if (
std::fabs(coeff[1] * coeff[5] - coeff[2] * coeff[4] - 1.)
723 double x0, y0, x1, y1;
727 for (
size_t i = 0; i < sourceSubCat.
size(); i++) {
728 x0 = sourceSubCat[i].getX();
729 y0 = sourceSubCat[i].getY();
731 auto refObjDist = searchNearestPoint(posRefSubCat, x1, y1,
733 if (refObjDist.first != posRefSubCat.
end()) {
738 std::cout <<
"Match: " << x0 <<
"," << y0 <<
" --> " 739 << x1 <<
"," << y1 <<
" <==> " 740 << refObjDist.first->getX() <<
"," << refObjDist.first->getY()
752 coeff = polyfit(1, srcMat, catMat);
755 for (
size_t i = 0; i < 6; ++i) {
761 matPair = FinalVerify(coeff, posRefProxyCat, sourceProxyCat,
764 std::cout <<
"Number of matches: " << matPair.size() <<
" vs " <<
771 if (matPair.size() > matPairSave.
size()) {
772 matPairSave = matPair;
780 if (matPairCand.
size() == 3) {
791 if (matPairCand.
size() == 0) {
794 size_t nmatch = matPairCand[0].
size();
796 for (
size_t i = 1; i < matPairCand.
size(); i++) {
797 if (matPairCand[i].size() > nmatch) {
798 nmatch = matPairCand[i].
size();
799 matPairRet = matPairCand[i];
std::int64_t RecordId
Type used for unique IDs for records.
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Eigen::Matrix2d getCdMatrix(Point2D const &pixel) const
Get the 2x2 CD matrix at the specified pixel position.
Extent< double, 2 > Extent2D
int numPointsForShape
"number of points in a matching shape" ;
AngleUnit constexpr radians
constant with units of radians
double matchingAllowancePix
"maximum allowed distance between reference objects and sources (pixels)" ;
std::vector< ReferenceMatch > ReferenceMatchVector
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Schema getSchema() const
Return the schema associated with the catalog's table.
ProxyVector makeProxies(lsst::afw::table::SourceCatalog const &sourceCat, afw::geom::SkyWcs const &distortedWcs, afw::geom::SkyWcs const &tanWcs)
std::vector< RecordProxy > ProxyVector
Include files required for standard LSST Exception handling.
A class representing an angle.
A base class for image defects.
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
iterator end()
Iterator access.
double maxDeterminant
"?" ;
Lightweight representation of a geometric match between two records.
const char * source()
Source function that allows astChannel to source from a Stream.
coord::IcrsCoord pixelToSky(Point2D const &pixel) const
Compute sky position(s) from pixel position(s)
Point2D skyToPixel(coord::IcrsCoord const &sky) const
Compute pixel position(s) from sky position(s)
table::Key< table::Array< double > > coeff
bool empty() const
Return true if the catalog has no records.
double maxOffsetPix
"maximum allowed frame translation (pixels)" ;
double degToRad(long double angleInDegrees)
Extent< double, 3 > Extent3D
#define LSST_EXCEPT(type,...)
Create an exception with a given type and message and optionally other arguments (dependent on the ty...
int numBrightStars
"maximum number of bright reference stars to use" ;
std::string refFluxField
"name of flux field in reference catalog" ;
Base::const_iterator const_iterator
A wrapper around a SimpleRecord or SourceRecord that allows us to record a pixel position in a way th...
lsst::afw::table::ReferenceMatchVector matchOptimisticB(lsst::afw::table::SimpleCatalog const &posRefCat, lsst::afw::table::SourceCatalog const &sourceCat, MatchOptimisticBControl const &control, afw::geom::SkyWcs const &wcs, int posRefBegInd=0, bool verbose=false)
Match sources to stars in a position reference catalog using optimistic pattern matching B...
double maxRotationDeg
"maximum allowed frame rotation (deg)" ;
size_type size() const
Return the number of elements in the catalog.
std::string sourceFluxField
"name of flux field in source catalog" ;
iterator begin()
Iterator access.
constexpr double asDegrees() const noexcept
Return an Angle's value in degrees.
table::Key< int > transform
double allowedNonperpDeg
"allowed non-perpendicularity of x and y axes (deg)" ;
A class to handle Icrs coordinates (inherits from Coord)
int minMatchedPairs
"minimum number of matches" ;
boost::shared_ptr< lsst::afw::table::SimpleRecord > record
A FunctorKey used to get or set celestial coordinates from a pair of Angle keys.