35 #include "Eigen/Eigenvalues"
46 template <
typename ImageT>
51 _constantWeight(constantWeight),
52 _eigenValues(
std::vector<double>()),
55 template <
typename ImageT>
57 template <
typename ImageT>
59 template <
typename ImageT>
61 template <
typename ImageT>
64 template <
typename ImageT>
66 if (_imageList.empty()) {
67 _dimensions = img->getDimensions();
69 if (getDimensions() != img->getDimensions()) {
72 img->getWidth() % img->getHeight() % _dimensions.getX() % _dimensions.getY())
81 _imageList.push_back(img);
85 template <
typename ImageT>
90 template <
typename ImageT>
92 if (_imageList.empty()) {
99 for (
typename ImageList::const_iterator
ptr = _imageList.begin(),
end = _imageList.end();
ptr !=
end;
103 *mean /= _imageList.size();
114 template <
typename T>
115 struct SortEvalueDecreasing
118 return a.first >
b.first;
123 template <
typename ImageT>
125 int const nImage = _imageList.size();
134 _eigenImages.clear();
137 _eigenValues.clear();
138 _eigenValues.push_back(1.0);
145 Eigen::MatrixXd R(nImage, nImage);
148 for (
int i = 0; i != nImage; ++i) {
150 double const flux_i = getFlux(i);
153 for (
int j = i; j != nImage; ++j) {
155 double const flux_j = getFlux(j);
158 if (_constantWeight) {
159 dot /= flux_i * flux_j;
161 R(i, j) = R(j, i) =
dot / nImage;
165 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(R);
166 Eigen::MatrixXd
const& Q = eVecValues.eigenvectors();
167 Eigen::VectorXd
const& lambda = eVecValues.eigenvalues();
173 lambdaAndIndex.
reserve(nImage);
175 for (
int i = 0; i != nImage; ++i) {
178 std::sort(lambdaAndIndex.
begin(), lambdaAndIndex.
end(), SortEvalueDecreasing<double>());
182 _eigenValues.clear();
183 _eigenValues.reserve(nImage);
184 for (
int i = 0; i != nImage; ++i) {
185 _eigenValues.push_back(lambdaAndIndex[i].
first);
192 _eigenImages.clear();
193 _eigenImages.reserve(ncomp < nImage ? ncomp : nImage);
195 for (
int i = 0; i < ncomp; ++i) {
200 int const ii = lambdaAndIndex[i].second;
205 for (
int j = 0; j != nImage; ++j) {
206 int const jj = lambdaAndIndex[j].second;
207 double const weight = Q(jj, ii) * (_constantWeight ? flux_bar / getFlux(jj) : 1);
208 eImage->scaledPlus(
weight, *_imageList[jj]);
210 _eigenImages.push_back(eImage);
220 template <
typename MaskedImageT>
224 MaskedImageT
const&
image
226 typedef typename MaskedImageT::Image ImageT;
230 }
else if (nEigen >
static_cast<int>(eigenImages.
size())) {
232 (
boost::format(
"You only have %d eigen images (you asked for %d)") %
233 eigenImages.
size() % nEigen)
240 Eigen::MatrixXd A(nEigen, nEigen);
241 Eigen::VectorXd
b(nEigen);
243 for (
int i = 0; i != nEigen; ++i) {
246 for (
int j = i; j != nEigen; ++j) {
247 A(i, j) = A(j, i) =
innerProduct(*eigenImages[i]->getImage(), *eigenImages[j]->getImage());
250 Eigen::VectorXd
x(nEigen);
253 x(0) =
b(0) / A(0, 0);
255 x = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(
b);
262 for (
int i = 0; i != nEigen; ++i) {
263 bestFitImage->scaledPlus(
x[i], *eigenImages[i]->getImage());
269 template <
typename ImageT>
272 unsigned long,
int const) {
276 template <
typename ImageT>
277 double do_updateBadPixels(detail::MaskedImage_tag
const&,
284 int const nImage = imageList.size();
287 "Please provide at least one Image for me to update");
292 double maxChange = 0.0;
298 for (
int i = 0; i != nImage; ++i) {
299 double const flux_i = fluxes[i];
301 for (
int y = 0;
y != height; ++
y) {
302 typename ImageT::const_x_iterator iptr = imageList[i]->row_begin(
y);
304 for (
typename ImageT::Image::x_iterator mptr = mean.row_begin(
y),
end = mean.row_end(
y);
305 mptr !=
end; ++mptr, ++iptr, ++wptr) {
306 if (!(iptr.mask() &
mask) && iptr.variance() > 0.0) {
308 float const var = iptr.variance() / (flux_i * flux_i);
309 float const ivar = 1.0 / var;
312 *mptr += value * ivar;
322 for (
int y = 0;
y != height; ++
y) {
324 for (
typename ImageT::Image::x_iterator mptr = mean.row_begin(
y),
end = mean.row_end(
y);
325 mptr !=
end; ++mptr, ++wptr) {
336 for (
int i = 0; i != nImage; ++i) {
337 double const flux_i = fluxes[i];
339 for (
int y = 0;
y != height; ++
y) {
340 typename ImageT::x_iterator iptr = imageList[i]->row_begin(
y);
341 for (
typename ImageT::Image::x_iterator mptr = mean.row_begin(
y),
end = mean.row_end(
y);
342 mptr !=
end; ++mptr, ++iptr) {
343 if ((iptr.mask() &
mask)) {
344 double const delta = ::fabs(flux_i * (*mptr) - iptr.image());
345 if (delta > maxChange) {
348 iptr.image() = flux_i * (*mptr);
354 if (ncomp >
static_cast<int>(eigenImages.size())) {
356 (
boost::format(
"You only have %d eigen images (you asked for %d)") %
357 eigenImages.size() % ncomp)
361 for (
int i = 0; i != nImage; ++i) {
363 fitEigenImagesToImage(eigenImages, ncomp, *imageList[i]);
365 for (
int y = 0;
y != height; ++
y) {
366 typename ImageT::x_iterator iptr = imageList[i]->row_begin(
y);
367 for (
typename ImageT::Image::const_x_iterator fptr = fitted->row_begin(
y),
368 end = fitted->row_end(
y);
369 fptr !=
end; ++fptr, ++iptr) {
370 if (iptr.mask() &
mask) {
371 double const delta =
fabs(
static_cast<double>(*fptr) - iptr.image());
372 if (delta > maxChange) {
376 iptr.image() = *fptr;
386 template <
typename ImageT>
388 return do_updateBadPixels<ImageT>(
typename ImageT::image_category(), _imageList, _fluxList, _eigenImages,
393 template <
typename T,
typename U>
395 IsSame(T
const&, U
const&) {}
396 bool operator()() {
return false; }
399 template <
typename T>
400 struct IsSame<T, T> {
401 IsSame(T
const& im1, T
const& im2) : _same(im1.row_begin(0) == im2.row_begin(0)) {}
402 bool operator()() {
return _same; }
409 template <
typename Image1T,
typename Image2T>
410 bool imagesAreIdentical(Image1T
const& im1, Image2T
const& im2) {
411 return IsSame<Image1T, Image2T>(im1, im2)();
414 template <
typename Image1T,
typename Image2T>
415 double innerProduct(Image1T
const& lhs, Image2T
const& rhs,
int border) {
416 if (lhs.getWidth() <= 2 * border || lhs.getHeight() <= 2 * border) {
418 (
boost::format(
"All image pixels are in the border of width %d: %dx%d") % border %
419 lhs.getWidth() % lhs.getHeight())
427 if (imagesAreIdentical(lhs, rhs)) {
428 for (
int y = border;
y != lhs.getHeight() - border; ++
y) {
429 for (
typename Image1T::const_x_iterator lptr = lhs.row_begin(
y) + border,
430 lend = lhs.row_end(
y) - border;
431 lptr != lend; ++lptr) {
439 if (lhs.getDimensions() != rhs.getDimensions()) {
441 (
boost::format(
"Dimension mismatch: %dx%d v. %dx%d") % lhs.getWidth() %
442 lhs.getHeight() % rhs.getWidth() % rhs.getHeight())
446 for (
int y = border;
y != lhs.getHeight() - border; ++
y) {
447 typename Image2T::const_x_iterator rptr = rhs.row_begin(
y) + border;
448 for (
typename Image1T::const_x_iterator lptr = lhs.row_begin(
y) + border,
449 lend = lhs.row_end(
y) - border;
450 lptr != lend; ++lptr, ++rptr) {
451 double const tmp = (*lptr) * (*rptr);
466 #define INSTANTIATE(T) \
467 template class ImagePca<Image<T> >; \
468 template double innerProduct(Image<T> const&, Image<T> const&, int); \
469 template class ImagePca<MaskedImage<T> >;
471 #define INSTANTIATE2(T, U) \
472 template double innerProduct(Image<T> const&, Image<U> const&, int); \
473 template double innerProduct(Image<U> const&, Image<T> const&, int);