Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DM-41634: Add PixelMap subclass using AST mappings #15

Merged
merged 6 commits into from
Dec 4, 2023
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
12 changes: 12 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,18 @@ add_compile_definitions(USE_YAML=TRUE)

find_library(FFTW3 fftw3 REQUIRED)

# find astshim library using EUPS package environent variable as hint
find_library(ASTSHIM_LIBRARIES
astshim REQUIRED
HINTS $ENV{ASTSHIM_DIR}/lib
)
#find astshim include directory
find_path(
ASTSHIM_INCLUDES
astshim.h REQUIRED
HINTS $ENV{ASTSHIM_DIR}/include
)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can't comment on this line below
This adds the target WCSfit as an executable
linking all static libs. astrometry which links astshim
Linking everything into one shlib gbdes might work and link this to this executable. Or remove
this target if it is not needed.

Suggested change
# add_subdirectory(src)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we need an executable target here, unless @cmsaunders still finds it convenient for comparison-testing sometimes.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm going to just leave in instructions for how to build the executable if needed.

# TMV or Eigen must be available, but TMV is preferred
find_package(TMV QUIET)
if(TMV_FOUND)
Expand Down
2 changes: 1 addition & 1 deletion astrometry/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
add_library(astrometry src/Astrometry.cpp src/PiecewiseMap.cpp src/PixelMapCollection.cpp src/SubMap.cpp
src/Wcs.cpp src/PixelMap.cpp src/PolyMap.cpp src/TemplateMap.cpp src/YAMLCollector.cpp)

target_include_directories(astrometry PUBLIC ${PROJECT_SOURCE_DIR}/astrometry/include)
target_include_directories(astrometry PUBLIC ${PROJECT_SOURCE_DIR}/astrometry/include ${ASTSHIM_INCLUDES})
target_link_libraries(astrometry PUBLIC Eigen3::Eigen ${FFTW3_LIBRARIES} yaml-cpp gbutil)
34 changes: 34 additions & 0 deletions astrometry/include/PixelMap.h
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,8 @@
#ifdef USE_YAML
#include "yaml-cpp/yaml.h"
#endif
#include "astshim.h"

#include "Std.h"
Copy link

@mwittgen mwittgen Nov 28, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is ugly in the original gbdes code that Std.h
Include all std headers and adding a using std

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, it's not great, but we decided not to go through and change all the original code. New code is supposed to ignore this, but I forgot in the instance you commented on below.

#include "LinearAlgebra.h"
#include "Astrometry.h"
Expand Down Expand Up @@ -332,5 +334,37 @@ namespace astrometry {
void operator=(const ColorTerm& rhs) = delete;

};

class ASTMap: public PixelMap {
// Class representing a constant shift of position
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand why the ASTMap class is described as representing solely a constant shift.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a mistake due to my using a different class as a template.

public:
ASTMap(const ast::Mapping& mapping_, string name_=""):

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is std::string obfuscated by the using std

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for catching this. We've been trying to include std:: in all new code, but I forgot and copied what was elsewhere in this file.

PixelMap(name_), mapping(mapping_) {}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be much better to make the mapping data member a std::shared_ptr<ast::Mapping> and call copy() here. You don't actually care about being able to share pointers, and all else being equal I'd prefer unique_ptr, but once you have a shared_ptr (and that's what methods like Mapping::copy return) you're stuck with it. The important thing is that the mapping_ could be an instance of a derived class of Mapping, and by copy-constructing it you're slicing off any extra stuff the derived class might have added.

Normally that would result in really horrible bugs, and the copy constructor would have been disabled to protect against that, and I find it strange that that's not the case in Mapping. I'm sure it's related to the fact that the real inheritance is happening in AST C code and the C++ inheritance is just a facade around that, but I still think it's much safer to hold by pointer and avoid the type-slicing yourself.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I've implemented what you suggest. Can you check and see if I did what you meant? I also changed duplicate so that it no longer does a copy.


virtual PixelMap* duplicate() const;

~ASTMap() {}

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can be declared as default dtor

Suggested change
~ASTMap() {}
~ASTMap() = default;


static string type() {return "AST";}

ast::Mapping mapping;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be better for this to be private.


void toWorld(double xpix, double ypix,
double &xworld, double &yworld,
double color=astrometry::NODATA) const;
void toPix(double xworld, double yworld,
double &xpix, double &ypix,
double color = astrometry::NODATA) const;

#ifdef USE_YAML
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you know if USE_YAML is actually set in our builds? If not, it might be cleaner to leaf this block out entirely.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, yaml is built into the mappings and if USE_YAML is false some key classes don't get built. Long term, it would be nice to remove this dependence though.

// AST maps are not expected to be used for output, so we do not
// attempt to write a full description to yaml.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could, because AST mappings can be serialized to ASDF YAML.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good to know. I'm not going to implement that now, because I think it's pretty unlikely it will be needed.

virtual void write(YAML::Emitter &os) const
{
os << YAML::BeginMap << YAML::Key << "Type" << YAML::Value << type()
<< YAML::EndMap;
}
#endif
};
} // namespace astrometry
#endif //PIXMAP_H
27 changes: 27 additions & 0 deletions astrometry/src/PixelMap.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -350,3 +350,30 @@ ConstantMap::create(const YAML::Node& node,
}
}
#endif

PixelMap*
ASTMap::duplicate() const {
std::shared_ptr<ast::Mapping> newMapping = mapping.copy();
return new ASTMap(*newMapping);
}

void
ASTMap::toWorld(double xpix, double ypix,
double& xworld, double& yworld,
double color) const {

std::vector<double> pixCoords{xpix, ypix};
std::vector<double> worldCoord = mapping.applyForward(pixCoords);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm a bit worried that this is going to be slow. AST is much faster if you give it many coordinates to transform and it's even faster if you ask it to transform a regular pixel grid to WCS in one go (astMapGrid).

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did a timing test, and there is a slow-down, but it's not prohibitive. This will only be used in a couple places in setting up the astrometric fit. For 31 visits, the time it took to do source association went from ~4 seconds to ~7 seconds, and the fit initialization went from ~33 seconds to ~41 seconds.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok. I assume it's not possible to run this in bulk given the way gbdes calls the transformation code, but I imagine it would pull all that time back if you did that.

xworld = worldCoord[0];
yworld = worldCoord[1];
}

void
ASTMap::toPix( double xworld, double yworld,
double &xpix, double &ypix,
double color) const {
std::vector<double> worldCoords(xworld, yworld);
std::vector<double> pixCoord = mapping.applyInverse(worldCoords);
xpix = pixCoord[0];
ypix = pixCoord[1];
}
4 changes: 2 additions & 2 deletions pydir/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@ set(PYTHON_MODULE_EXTENSION ".so")

pybind11_add_module(wcsfit wcsfit.cpp)

target_link_libraries(wcsfit PUBLIC subs cfitsio)
target_link_libraries(wcsfit PUBLIC subs cfitsio ${ASTSHIM_LIBRARIES})

target_include_directories(wcsfit PUBLIC ${PROJECT_SOURCE_DIR}/include)
target_include_directories(wcsfit PUBLIC ${PROJECT_SOURCE_DIR}/include ${ASTSHIM_INCLUDES})

install(TARGETS wcsfit DESTINATION "${PROJECT_SOURCE_DIR}")
5 changes: 5 additions & 0 deletions pydir/wcsfit.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include "TPVMap.h"
#include "WCSFoFRoutine.h"
#include "WCSFitRoutine.h"
#include "astshim.h"

using namespace pybind11::literals;
namespace py = pybind11;
Expand Down Expand Up @@ -102,6 +103,10 @@ PYBIND11_MODULE(wcsfit, m) {
py::class_<astrometry::SubMap, astrometry::PixelMap>(m, "SubMap")
.def(py::init<list<astrometry::PixelMap *> const &, std::string, bool>());

py::class_<astrometry::ASTMap, astrometry::PixelMap>(m, "ASTMap")
.def(py::init<ast::Mapping, std::string>(),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
.def(py::init<ast::Mapping, std::string>(),
.def(py::init<ast::Mapping const &, std::string>(),

If this isn't marked as a reference or pointer we'll get more type-slicing in the pybind11 layer.

py::arg("mapping_"), py::arg("name_") = "");

py::class_<astrometry::Wcs, std::shared_ptr<astrometry::Wcs>>(m, "Wcs")
.def(py::init<astrometry::PixelMap *, astrometry::SphericalCoords const &, std::string, double,
bool>(),
Expand Down
1 change: 1 addition & 0 deletions ups/gbdes.table
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
setupRequired(astshim)
envPrepend(PYTHONPATH, ${PRODUCT_DIR}/build/pydir)