34 #include "Eigen/Eigenvalues"
45 template <
typename ImageT>
50 _constantWeight(constantWeight),
51 _eigenValues(
std::vector<double>()),
54 template <
typename ImageT>
56 template <
typename ImageT>
58 template <
typename ImageT>
60 template <
typename ImageT>
63 template <
typename ImageT>
65 if (_imageList.empty()) {
66 _dimensions = img->getDimensions();
68 if (getDimensions() != img->getDimensions()) {
70 (
boost::format(
"Dimension mismatch: saw %dx%d; expected %dx%d") %
71 img->getWidth() % img->getHeight() % _dimensions.getX() % _dimensions.getY())
80 _imageList.push_back(img);
81 _fluxList.push_back(flux);
84 template <
typename ImageT>
89 template <
typename ImageT>
91 if (_imageList.empty()) {
98 for (
typename ImageList::const_iterator
ptr = _imageList.begin(),
end = _imageList.end();
ptr !=
end;
102 *mean /= _imageList.size();
113 template <
typename T>
114 struct SortEvalueDecreasing
117 return a.first >
b.first;
122 template <
typename ImageT>
124 int const nImage = _imageList.size();
133 _eigenImages.clear();
136 _eigenValues.clear();
137 _eigenValues.push_back(1.0);
144 Eigen::MatrixXd R(nImage, nImage);
147 for (
int i = 0; i != nImage; ++i) {
149 double const flux_i = getFlux(i);
152 for (
int j = i; j != nImage; ++j) {
154 double const flux_j = getFlux(j);
157 if (_constantWeight) {
158 dot /= flux_i * flux_j;
160 R(i, j) = R(j, i) =
dot / nImage;
164 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(R);
165 Eigen::MatrixXd
const& Q = eVecValues.eigenvectors();
166 Eigen::VectorXd
const& lambda = eVecValues.eigenvalues();
172 lambdaAndIndex.
reserve(nImage);
174 for (
int i = 0; i != nImage; ++i) {
177 std::sort(lambdaAndIndex.
begin(), lambdaAndIndex.
end(), SortEvalueDecreasing<double>());
181 _eigenValues.clear();
182 _eigenValues.reserve(nImage);
183 for (
int i = 0; i != nImage; ++i) {
184 _eigenValues.push_back(lambdaAndIndex[i].
first);
191 _eigenImages.clear();
192 _eigenImages.reserve(ncomp < nImage ? ncomp : nImage);
194 for (
int i = 0; i < ncomp; ++i) {
199 int const ii = lambdaAndIndex[i].second;
204 for (
int j = 0; j != nImage; ++j) {
205 int const jj = lambdaAndIndex[j].second;
206 double const weight = Q(jj, ii) * (_constantWeight ? flux_bar / getFlux(jj) : 1);
207 eImage->scaledPlus(
weight, *_imageList[jj]);
209 _eigenImages.push_back(eImage);
219 template <
typename MaskedImageT>
223 MaskedImageT
const&
image
225 using ImageT =
typename MaskedImageT::Image;
229 }
else if (nEigen >
static_cast<int>(eigenImages.
size())) {
231 (
boost::format(
"You only have %d eigen images (you asked for %d)") %
232 eigenImages.
size() % nEigen)
239 Eigen::MatrixXd A(nEigen, nEigen);
240 Eigen::VectorXd
b(nEigen);
242 for (
int i = 0; i != nEigen; ++i) {
245 for (
int j = i; j != nEigen; ++j) {
246 A(i, j) = A(j, i) =
innerProduct(*eigenImages[i]->getImage(), *eigenImages[j]->getImage());
249 Eigen::VectorXd
x(nEigen);
252 x(0) =
b(0) / A(0, 0);
254 x = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(
b);
261 for (
int i = 0; i != nEigen; ++i) {
262 bestFitImage->scaledPlus(
x[i], *eigenImages[i]->getImage());
268 template <
typename ImageT>
271 unsigned long,
int const) {
275 template <
typename ImageT>
276 double do_updateBadPixels(detail::MaskedImage_tag
const&,
283 int const nImage = imageList.size();
286 "Please provide at least one Image for me to update");
291 double maxChange = 0.0;
297 for (
int i = 0; i != nImage; ++i) {
298 double const flux_i = fluxes[i];
300 for (
int y = 0;
y != height; ++
y) {
301 typename ImageT::const_x_iterator iptr = imageList[i]->row_begin(
y);
303 for (
typename ImageT::Image::x_iterator mptr = mean.row_begin(
y),
end = mean.row_end(
y);
304 mptr !=
end; ++mptr, ++iptr, ++wptr) {
305 if (!(iptr.mask() &
mask) && iptr.variance() > 0.0) {
307 float const var = iptr.variance() / (flux_i * flux_i);
308 float const ivar = 1.0 / var;
311 *mptr += value * ivar;
321 for (
int y = 0;
y != height; ++
y) {
323 for (
typename ImageT::Image::x_iterator mptr = mean.row_begin(
y),
end = mean.row_end(
y);
324 mptr !=
end; ++mptr, ++wptr) {
335 for (
int i = 0; i != nImage; ++i) {
336 double const flux_i = fluxes[i];
338 for (
int y = 0;
y != height; ++
y) {
339 typename ImageT::x_iterator iptr = imageList[i]->row_begin(
y);
340 for (
typename ImageT::Image::x_iterator mptr = mean.row_begin(
y),
end = mean.row_end(
y);
341 mptr !=
end; ++mptr, ++iptr) {
342 if ((iptr.mask() &
mask)) {
343 double const delta = ::fabs(flux_i * (*mptr) - iptr.image());
344 if (delta > maxChange) {
347 iptr.image() = flux_i * (*mptr);
353 if (ncomp >
static_cast<int>(eigenImages.size())) {
355 (
boost::format(
"You only have %d eigen images (you asked for %d)") %
356 eigenImages.size() % ncomp)
360 for (
int i = 0; i != nImage; ++i) {
362 fitEigenImagesToImage(eigenImages, ncomp, *imageList[i]);
364 for (
int y = 0;
y != height; ++
y) {
365 typename ImageT::x_iterator iptr = imageList[i]->row_begin(
y);
366 for (
typename ImageT::Image::const_x_iterator fptr = fitted->row_begin(
y),
367 end = fitted->row_end(
y);
368 fptr !=
end; ++fptr, ++iptr) {
369 if (iptr.mask() &
mask) {
370 double const delta =
fabs(
static_cast<double>(*fptr) - iptr.image());
371 if (delta > maxChange) {
375 iptr.image() = *fptr;
385 template <
typename ImageT>
387 return do_updateBadPixels<ImageT>(
typename ImageT::image_category(), _imageList, _fluxList, _eigenImages,
392 template <
typename T,
typename U>
394 IsSame(
T const&, U
const&) {}
395 bool operator()() {
return false; }
398 template <
typename T>
399 struct IsSame<
T,
T> {
400 IsSame(
T const& im1,
T const& im2) : _same(im1.row_begin(0) == im2.row_begin(0)) {}
401 bool operator()() {
return _same; }
408 template <
typename Image1T,
typename Image2T>
409 bool imagesAreIdentical(Image1T
const& im1, Image2T
const& im2) {
410 return IsSame<Image1T, Image2T>(im1, im2)();
413 template <
typename Image1T,
typename Image2T>
414 double innerProduct(Image1T
const& lhs, Image2T
const& rhs,
int border) {
415 if (lhs.getWidth() <= 2 * border || lhs.getHeight() <= 2 * border) {
417 (
boost::format(
"All image pixels are in the border of width %d: %dx%d") % border %
418 lhs.getWidth() % lhs.getHeight())
426 if (imagesAreIdentical(lhs, rhs)) {
427 for (
int y = border;
y != lhs.getHeight() - border; ++
y) {
428 for (
typename Image1T::const_x_iterator lptr = lhs.row_begin(
y) + border,
429 lend = lhs.row_end(
y) - border;
430 lptr != lend; ++lptr) {
438 if (lhs.getDimensions() != rhs.getDimensions()) {
440 (
boost::format(
"Dimension mismatch: %dx%d v. %dx%d") % lhs.getWidth() %
441 lhs.getHeight() % rhs.getWidth() % rhs.getHeight())
445 for (
int y = border;
y != lhs.getHeight() - border; ++
y) {
446 typename Image2T::const_x_iterator rptr = rhs.row_begin(
y) + border;
447 for (
typename Image1T::const_x_iterator lptr = lhs.row_begin(
y) + border,
448 lend = lhs.row_end(
y) - border;
449 lptr != lend; ++lptr, ++rptr) {
450 double const tmp = (*lptr) * (*rptr);
465 #define INSTANTIATE(T) \
466 template class ImagePca<Image<T> >; \
467 template double innerProduct(Image<T> const&, Image<T> const&, int); \
468 template class ImagePca<MaskedImage<T> >;
470 #define INSTANTIATE2(T, U) \
471 template double innerProduct(Image<T> const&, Image<U> const&, int); \
472 template double innerProduct(Image<U> const&, Image<T> const&, int);
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
afw::table::PointKey< int > dimensions
#define INSTANTIATE2(ImagePixelT1, ImagePixelT2)
A class to represent a 2-dimensional array of pixels.
virtual double updateBadPixels(unsigned long mask, int const ncomp)
Update the bad pixels (i.e.
void addImage(std::shared_ptr< ImageT > img, double flux=0.0)
Add an image to the set to be analyzed.
ImageList getImageList() const
Return the list of images being analyzed.
ImagePca(bool constantWeight=true)
ctor
std::shared_ptr< ImageT > getMean() const
Return the mean of the images in ImagePca's list.
std::vector< std::shared_ptr< ImageT > > ImageList
ImagePca & operator=(ImagePca const &)
Reports attempts to exceed implementation-defined length limits for some classes.
Reports attempts to access elements outside a valid range of indices.
T emplace_back(T... args)
#define INSTANTIATE(TYPE)
def dot(symb, c, r, frame=None, size=2, ctype=None, origin=afwImage.PARENT, *args, **kwargs)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
double innerProduct(Image1T const &lhs, Image2T const &rhs, int const border=0)
Calculate the inner product of two images.
float Pixel
Typedefs to be used for pixel values.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
A base class for image defects.
afw::table::Key< double > weight