diff --git a/include/IRFFT.h b/include/IRFFT.h new file mode 100644 index 0000000..d028ab1 --- /dev/null +++ b/include/IRFFT.h @@ -0,0 +1,407 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : fft.hxx +// Author : Pavel A. Koshevoy +// Created : 2005/06/03 10:16 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : Wrapper class and helper functions for working with +// FFTW3 and ITK images. + +#ifndef FFT_HXX_ +#define FFT_HXX_ + +// FFTW includes: +#include + +// ITK includes: +#include + +// system includes: +#include +#include +#include + + +namespace itk_fft +{ + //---------------------------------------------------------------- + // set_num_fftw_threads + // + // set's number of threads used by fftw, returns previous value: + extern std::size_t set_num_fftw_threads(std::size_t num_threads); + + //---------------------------------------------------------------- + // itk_image_t + // + typedef itk::Image itk_image_t; + + //---------------------------------------------------------------- + // itk_imageptr_t + // + typedef itk_image_t::Pointer itk_imageptr_t; + + //---------------------------------------------------------------- + // fft_complex_t + // + typedef std::complex fft_complex_t; + + //---------------------------------------------------------------- + // fft_data_t + // + class fft_data_t + { + public: + fft_data_t(): data_(NULL), nx_(0), ny_(0) {} + fft_data_t(const unsigned int w, const unsigned int h); + explicit fft_data_t(const itk_imageptr_t & real); + fft_data_t(const itk_imageptr_t & real, + const itk_imageptr_t & imag); + fft_data_t(const fft_data_t & data); + ~fft_data_t() { cleanup(); } + + fft_data_t & operator = (const fft_data_t & data); + + void cleanup(); + + void resize(const unsigned int w, const unsigned int h); + + void fill(const float real, const float imag = 0.0); + + void setup(const itk_imageptr_t & real, + const itk_imageptr_t & imag = itk_imageptr_t(NULL)); + + // ITK helpers: + itk_imageptr_t component(const bool imag = 0) const; + inline itk_imageptr_t real() const { return component(0); } + inline itk_imageptr_t imag() const { return component(1); } + + // Apply a low-pass filter to this image. This function will + // zero-out high-frequency components, where the cutoff frequency + // is specified by radius r in [0, 1]. The sharpness of the + // cutoff may be controlled by parameter s, where s == 0 results in + // an ideal low-pass filter, and s == 1 is a low pass filter defined + // by a scaled and shifted cosine function, 1 at the origin, + // 0.5 at the cutoff frequency and 0 at twice the cutoff frequency. + void apply_lp_filter(const double r, const double s = 0); + + // accessors: + inline const fft_complex_t * data() const { return data_; } + inline fft_complex_t * data() { return data_; } + + inline unsigned int nx() const { return nx_; } + inline unsigned int ny() const { return ny_; } + + inline const fft_complex_t & operator() (const unsigned int & x, + const unsigned int & y) const + { return at(x, y); } + + inline fft_complex_t & operator() (const unsigned int & x, + const unsigned int & y) + { return at(x, y); } + + inline const fft_complex_t & at(const unsigned int & x, + const unsigned int & y) const + { return data_[y + ny_ * x]; } + + inline fft_complex_t & at(const unsigned int & x, + const unsigned int & y) + { return data_[y + ny_ * x]; } + + protected: + fft_complex_t * data_; + unsigned int nx_; + unsigned int ny_; + }; + + + //---------------------------------------------------------------- + // fft + // + extern bool + fft(itk_image_t::ConstPointer & in, fft_data_t & out); + + //---------------------------------------------------------------- + // fft + // + extern bool + fft(const fft_data_t & in, fft_data_t & out, int sign = FFTW_FORWARD); + + //---------------------------------------------------------------- + // ifft + // + extern bool + ifft(const fft_data_t & in, fft_data_t & out); + + //---------------------------------------------------------------- + // ifft + // + inline fft_data_t + ifft(const fft_data_t & in) + { + fft_data_t out; +#ifndef NDEBUG // get around an annoying compiler warning: + bool ok = +#endif + ifft(in, out); + assert(ok); + return out; + } + + + //---------------------------------------------------------------- + // fn_fft_c_t + // + typedef fft_complex_t(*fn_fft_c_t)(const fft_complex_t & in); + + //---------------------------------------------------------------- + // fn_fft_cc_t + // + typedef fft_complex_t(*fn_fft_cc_t)(const fft_complex_t & a, + const fft_complex_t & b); + + //---------------------------------------------------------------- + // elem_by_elem + // + extern void + elem_by_elem(fn_fft_c_t f, + const fft_data_t & in, + fft_data_t & out); + + //---------------------------------------------------------------- + // elem_by_elem + // + extern void + elem_by_elem(fft_complex_t(*f)(const float & a, + const fft_complex_t & b), + const float & a, + const fft_data_t & b, + fft_data_t & c); + + //---------------------------------------------------------------- + // elem_by_elem + // + extern void + elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, + const float & b), + const fft_data_t & a, + const float & b, + fft_data_t & c); + + //---------------------------------------------------------------- + // elem_by_elem + // + extern void + elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, + const fft_complex_t & b), + const fft_complex_t & a, + const fft_data_t & b, + fft_data_t & c); + + //---------------------------------------------------------------- + // elem_by_elem + // + extern void + elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, + const fft_complex_t & b), + const fft_data_t & a, + const fft_complex_t & b, + fft_data_t & c); + + //---------------------------------------------------------------- + // elem_by_elem + // + extern void + elem_by_elem(fn_fft_cc_t f, + const fft_data_t & a, + const fft_data_t & b, + fft_data_t & c); + + + //---------------------------------------------------------------- + // conj + // + // element-by-element complex conjugate: + // + inline void + conj(const fft_data_t & in, fft_data_t & out) + { elem_by_elem(std::conj, in, out); } + + //---------------------------------------------------------------- + // conj + // + inline fft_data_t + conj(const fft_data_t & in) + { + fft_data_t out; + conj(in, out); + return out; + } + + + //---------------------------------------------------------------- + // sqrt + // + // element-by-element square root: + // + inline void + sqrt(const fft_data_t & in, fft_data_t & out) + { elem_by_elem(std::sqrt, in, out); } + + //---------------------------------------------------------------- + // sqrt + // + inline fft_data_t + sqrt(const fft_data_t & in) + { + fft_data_t out; + sqrt(in, out); + return out; + } + + + //---------------------------------------------------------------- + // _mul + // + template + inline fft_complex_t + _mul(const a_t & a, const b_t & b) + { return a * b; } + + //---------------------------------------------------------------- + // mul + // + // element-by-element multiplication, c = a * b: + // + template + inline void + mul(const a_t & a, const b_t & b, fft_data_t & c) + { elem_by_elem(_mul, a, b, c); } + + //---------------------------------------------------------------- + // mul + // + template + inline fft_data_t + mul(const a_t & a, const b_t & b) + { + fft_data_t c; + mul(a, b, c); + return c; + } + + + //---------------------------------------------------------------- + // _div + // + template + inline fft_complex_t + _div(const a_t & a, const b_t & b) + { return a / b; } + + //---------------------------------------------------------------- + // div + // + // element-by-element division, c = a / b: + // + template + inline void + div(const a_t & a, const b_t & b, fft_data_t & c) + { elem_by_elem(_div, a, b, c); } + + //---------------------------------------------------------------- + // div + // + template + inline fft_data_t + div(const a_t & a, const b_t & b) + { + fft_data_t c; + div(a, b, c); + return c; + } + + + //---------------------------------------------------------------- + // _add + // + template + inline fft_complex_t + _add(const a_t & a, const b_t & b) + { return a + b; } + + //---------------------------------------------------------------- + // add + // + // element-by-element addition, c = a + b: + // + template + inline void + add(const a_t & a, const b_t & b, fft_data_t & c) + { elem_by_elem(_add, a, b, c); } + + //---------------------------------------------------------------- + // add + // + template + inline fft_data_t + add(const a_t & a, const b_t & b) + { + fft_data_t c; + add(a, b, c); + return c; + } + + + //---------------------------------------------------------------- + // _sub + // + template + inline fft_complex_t + _sub(const a_t & a, const b_t & b) + { return a - b; } + + //---------------------------------------------------------------- + // sub + // + // element-by-element subtraction, c = a - b: + // + template + inline void + sub(const a_t & a, const b_t & b, fft_data_t & c) + { elem_by_elem(_sub, a, b, c); } + + //---------------------------------------------------------------- + // sub + // + template + inline fft_data_t + sub(const a_t & a, const b_t & b) + { + fft_data_t c; + sub(a, b, c); + return c; + } +} + + +#endif // FFT_HXX_ diff --git a/include/IRFFTCommon.h b/include/IRFFTCommon.h new file mode 100644 index 0000000..6c4db1c --- /dev/null +++ b/include/IRFFTCommon.h @@ -0,0 +1,594 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : fft_common.hxx +// Author : Pavel A. Koshevoy +// Created : 2005/11/10 14:05 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : Helper functions for image alignment (registration) +// using phase correlation to find the translation vector. + +#ifndef FFT_COMMON_HXX_ +#define FFT_COMMON_HXX_ + +// local includes: +#include "itkIRCommon.h" +#include "IRFFT.h" +#include "itkIRDynamicArray.h" +#include "IRLog.h" + +// system includes: +#include +#include +#include + +#ifndef WIN32 +# include +#endif + +// namespace access: +using namespace itk_fft; + +//---------------------------------------------------------------- +// DEBUG_PDF +// +// #define DEBUG_PDF + +#ifdef DEBUG_PDF +// FIXME: this is not thread safe: +# define DEBUG_PDF_PFX the_text_t("PDF-") +# define DEBUG_CLUSTERS +# define DEBUG_MARKERS +extern unsigned int DEBUG_COUNTER1; +extern unsigned int DEBUG_COUNTER2; +#else +# define DEBUG_PDF_PFX the_text_t("") +#endif + + +//---------------------------------------------------------------- +// local_max_t +// +class local_max_t +{ +public: + local_max_t(const double value, const double x, const double y, const unsigned int area) + : value_(value) + , x_(x) + , y_(y) + , area_(area) + {} + + inline bool + operator==(const local_max_t & lm) const + { + return (value_ == lm.value_) && (x_ == lm.x_) && (y_ == lm.y_); + } + + inline bool + operator<(const local_max_t & lm) const + { + return value_ < lm.value_; + } + + inline bool + operator>(const local_max_t & lm) const + { + return value_ > lm.value_; + } + + // cluster value: + double value_; + + // center of mass of the cluster: + double x_; + double y_; + + // area of the cluster: + unsigned int area_; +}; + +//---------------------------------------------------------------- +// operator << +// +inline std::ostream & +operator<<(std::ostream & sout, const local_max_t & lm) +{ + return sout << lm.value_ << '\t' << lm.x_ << '\t' << lm.y_ << '\t' << lm.area_ << endl; +} + +//---------------------------------------------------------------- +// cluster_bbox_t +// +class cluster_bbox_t +{ +public: + cluster_bbox_t() { reset(); } + + void + reset() + { + min_[0] = std::numeric_limits::max(); + min_[1] = std::numeric_limits::max(); + max_[0] = std::numeric_limits::min(); + max_[1] = std::numeric_limits::min(); + } + + void + update(int x, int y) + { + min_[0] = std::min(min_[0], x); + min_[1] = std::min(min_[1], y); + max_[0] = std::max(max_[0], x); + max_[1] = std::max(max_[1], y); + } + + int min_[2]; + int max_[2]; +}; + +//---------------------------------------------------------------- +// find_maxima_cm +// +// The percentage below refers to the number of pixels that fall +// below the maxima. Thus, the number of pixels above the maxima +// is 1 - percentage. This way it is possible to specify a +// thresholding value without knowing anything about the image. +// +// The given image is thresholded, the resulting clusters/blobs +// are identified/classified, the center of mass of each cluster +// is treated as a maxima in the image. +// +// Returns the number of maxima found. +// +extern unsigned int +find_maxima_cm(std::list & max_list, + const itk_image_t::Pointer & image, + const double percentage = 0.9995, + const the_text_t & prefix = DEBUG_PDF_PFX, + const the_text_t & suffix = the_text_t(".png")); + + +//---------------------------------------------------------------- +// find_correlation +// +template +unsigned int +find_correlation(std::list & max_list, + const TImage * fi, + const TImage * mi, + double lp_filter_r, + double lp_filter_s, + const double min_overlap, + const double max_overlap) +{ + itk_image_t::Pointer z0 = cast(fi); + itk_image_t::Pointer z1 = cast(mi); + return find_correlation(max_list, z0, z1, lp_filter_r, lp_filter_s, min_overlap, max_overlap); +} + + +//---------------------------------------------------------------- +// find_correlation +// +template <> +unsigned int +find_correlation(std::list & max_list, + const itk_image_t * fi, + const itk_image_t * mi, + + + // low pass filter parameters + // (resampled data requires less smoothing): + double lp_filter_r, + double lp_filter_s, + const double overlap_min, + const double overlap_max); + + +//---------------------------------------------------------------- +// find_correlation +// +// This is a backwards compatibility API +// +template +unsigned int +find_correlation(const TImage * fi, + const TImage * mi, + std::list & max_list, + bool resampled_data, + const double overlap_min, + const double overlap_max) +{ + double lp_filter_r = resampled_data ? 0.9 : 0.5; + return find_correlation(max_list, fi, mi, lp_filter_r, 0.1, overlap_min, overlap_max); +} + + +//---------------------------------------------------------------- +// threshold_maxima +// +// Discard maxima whose mass is below a given threshold ratio +// of the total mass of all maxima: +// +void +threshold_maxima(std::list & max_list, const double threshold); + +//---------------------------------------------------------------- +// reject_negligible_maxima +// +// Discard maxima that are worse than the best maxima by a factor +// greater than the given threshold ratio: +// +// Returns the size of the new list. +// +unsigned int +reject_negligible_maxima(std::list & max_list, const double threshold); + +//---------------------------------------------------------------- +// overlap_t +// +class overlap_t +{ +public: + overlap_t() + : id_(~0) + , overlap_(0.0) + {} + + overlap_t(const unsigned int id, const double overlap) + : id_(id) + , overlap_(overlap) + {} + + inline bool + operator==(const overlap_t & d) const + { + return id_ == d.id_; + } + + inline bool + operator<(const overlap_t & d) const + { + return overlap_ < d.overlap_; + } + + inline bool + operator>(const overlap_t & d) const + { + return overlap_ > d.overlap_; + } + + unsigned int id_; + double overlap_; +}; + +//---------------------------------------------------------------- +// reject_negligible_overlap +// +void +reject_negligible_overlap(std::list & ol, const double threshold); + + +//---------------------------------------------------------------- +// estimate_displacement +// +template +double +estimate_displacement(the_log_t & log, + const TImage * a, + const TImage * b, + const local_max_t & lm, + translate_transform_t::Pointer & transform, + image_t::PointType offset_min, + image_t::PointType offset_max, + const double overlap_min = 0.0, + const double overlap_max = 1.0, + const mask_t * mask_a = NULL, + const mask_t * mask_b = NULL) +{ + // FIXME: +#ifdef DEBUG_PDF + the_text_t fn_debug = the_text_t::number(DEBUG_COUNTER1, 3, '0') + the_text_t("-") + + the_text_t::number(DEBUG_COUNTER2, 1, '0') + the_text_t("-"); +#endif + + itk_image_t::SizeType max_sz = calc_padding(a, b); + const unsigned int & w = max_sz[0]; + const unsigned int & h = max_sz[1]; + + // evaluate 4 permutations of the maxima: + // typedef itk::NormalizedCorrelationImageToImageMetric + typedef itk::MeanSquaresImageToImageMetric metric_t; + typedef itk::LinearInterpolateImageFunction interpolator_t; + + typedef typename TImage::SpacingType spacing_t; + spacing_t spacing = a->GetSpacing(); + double sx = spacing[0]; + double sy = spacing[1]; + + // evaluate 4 permutation of the maxima: + // TODO: James A: Is this for the fft? We can run the FFT in such a way that we just know... + const double x = lm.x_; + const double y = lm.y_; + + const vec2d_t t[] = { + vec2d(sx * x, sy * y), vec2d(sx * (x - w), sy * y), vec2d(sx * x, sy * (y - h)), vec2d(sx * (x - w), sy * (y - h)) + }; + + double best_metric = std::numeric_limits::max(); + itk_image_t::SizeType max_sz_with_spacing = max_sz; + max_sz_with_spacing[0] = max_sz[0] * sx; + max_sz_with_spacing[1] = max_sz[1] * sy; + + + for (unsigned int i = 0; i < 4; i++) + { + double overlap = OverlapPercent(max_sz_with_spacing, t[i]); + + if (overlap < overlap_min || overlap > overlap_max) + continue; + + translate_transform_t::Pointer ti = translate_transform_t::New(); + ti->SetOffset(t[i]); + + // FIXME: +#ifdef DEBUG_PDF + the_text_t fn = fn_debug + the_text_t::number(i) + ".png"; + save_rgb(fn, a, b, ti, mask_a, mask_b); +#endif + + const double area_ratio = overlap_ratio(a, b, ti); + int old_precision = log.precision(); + log.precision(2); + log << i << ": " << setw(3) << std::fixed << area_ratio * 100.0 << "% of overlap, "; + log.precision(old_precision); + + if (area_ratio < overlap_min || area_ratio > overlap_max) + { + log << "skipping..." << endl; + continue; + } + + // double metric = eval_metric(ti, a, b); + double metric = my_metric(a, b, ti, mask_a, mask_b, overlap_min, overlap_max); + log << std::scientific << metric; + if (metric < best_metric) + { + transform = ti; + best_metric = metric; + log << " - better..." << endl; + +#ifdef DEBUG_PDF + save_rgb(fn_debug + the_text_t::number(i) + "-better.png", a, b, ti, mask_a, mask_b); +#endif + } + else + { + log << " - worse..." << endl; + } + } + + return best_metric; +} + + +//---------------------------------------------------------------- +// match_one_pair +// +// match two images, find the best matching transform and +// return the corresponding match metric value +// +template +double +match_one_pair(the_log_t & log, + const bool images_were_resampled, + const bool use_std_mask, + const TImage * fi, + const TImage * mi, + const TMask * fi_mask, + const TMask * mi_mask, + + const double overlap_min, + const double overlap_max, + + image_t::PointType offset_min, + image_t::PointType offset_max, + + translate_transform_t::Pointer & ti, + std::list & peaks, + unsigned int & num_peaks, + + // ideally this should be one, but radial distortion may + // generate several valid peaks (up to 4), so it may be + // necessary to consider more peaks for the unmatched images: + const unsigned int max_peaks) +{ +#ifdef DEBUG_PDF + DEBUG_COUNTER1++; +#endif + + ti = NULL; + + unsigned int total_peaks = 0; + if (use_std_mask) + { + // ugh, blank out 5 percent of the image at the bottom + // to cover up the image id: + typename TImage::SizeType fi_sz = fi->GetLargestPossibleRegion().GetSize(); + typename TImage::SizeType mi_sz = mi->GetLargestPossibleRegion().GetSize(); + + unsigned int fi_y = (unsigned int)(0.95 * fi_sz[1]); + unsigned int mi_y = (unsigned int)(0.95 * mi_sz[1]); + + typename TImage::Pointer fi_filled = cast(fi); + fill(fi_filled, 0, fi_y, fi_sz[0], fi_sz[1] - fi_y, 0); + + typename TImage::Pointer mi_filled = cast(mi); + fill(mi_filled, 0, mi_y, mi_sz[0], mi_sz[1] - mi_y, 0); + + total_peaks = + find_correlation(fi_filled, mi_filled, peaks, images_were_resampled, overlap_min, overlap_max); + } + else + { + total_peaks = find_correlation(fi, mi, peaks, images_were_resampled, overlap_min, overlap_max); + } + + num_peaks = reject_negligible_maxima(peaks, 3.0); + log << num_peaks << '/' << total_peaks << " eligible peaks, "; + + if (num_peaks > max_peaks) + { + log << "skipping..." << endl; + return std::numeric_limits::max(); + } + + // choose the best peak: + double best_metric = std::numeric_limits::max(); + +#ifdef DEBUG_PDF + DEBUG_COUNTER2 = 0; +#endif + for (std::list::iterator j = peaks.begin(); j != peaks.end(); ++j) + { + log << "evaluating permutations..." << endl; + const local_max_t & lm = *j; + + translate_transform_t::Pointer tmp = translate_transform_t::New(); + double metric = estimate_displacement( + log, fi, mi, lm, tmp, offset_min, offset_max, overlap_min, overlap_max, fi_mask, mi_mask); + if (metric < best_metric) + { + best_metric = metric; + ti = tmp; + } + +#ifdef DEBUG_PDF + DEBUG_COUNTER2++; +#endif + } + + return best_metric; +} + + +//---------------------------------------------------------------- +// match_one_pair +// +// match two images, find the best matching transform and +// return the corresponding match metric value +// +template +double +match_one_pair(the_log_t & log, + const bool images_were_resampled, + const bool use_std_mask, + const TImage * fi, + const TImage * mi, + const TMask * fi_mask, + const TMask * mi_mask, + + const double overlap_min, + const double overlap_max, + + image_t::PointType offset_min, + image_t::PointType offset_max, + + const unsigned int node_id, + translate_transform_t::Pointer & ti, + std::pair> & peaks, + + // ideally this should be one, but radial distortion may + // generate several valid peaks (up to 4), so it may be + // necessary to consider more peaks for the unmatched images: + const unsigned int max_peaks) +{ + unsigned int peak_list_size = 0; + std::list peak_list; + double m = match_one_pair(log, + images_were_resampled, + use_std_mask, + fi, + mi, + fi_mask, + mi_mask, + overlap_min, + overlap_max, + offset_min, + offset_max, + ti, + peak_list, + peak_list_size, + max_peaks); + + // this info will be used when trying to match the unmatched images: + if (peak_list_size != 0 && peak_list_size <= max_peaks && + (peaks.second.empty() || peaks.second.size() > peak_list_size)) + { + peaks.first = node_id; + peaks.second.clear(); + peaks.second.splice(peaks.second.end(), peak_list); + } + + return m; +} + +//---------------------------------------------------------------- +// match_one_pair +// +template +double +match_one_pair(the_log_t & log, + const bool images_were_resampled, + const bool use_std_mask, + const TImage * fi, + const TImage * mi, + const TMask * fi_mask, + const TMask * mi_mask, + const double overlap_min, + const double overlap_max, + image_t::PointType offset_min, + image_t::PointType offset_max, + translate_transform_t::Pointer & ti) +{ + std::list peaks; + unsigned int num_peaks = 0; + return match_one_pair(log, + images_were_resampled, + use_std_mask, + fi, + mi, + fi_mask, + mi_mask, + overlap_min, + overlap_max, + offset_min, + offset_max, + ti, + peaks, + num_peaks, + UINT_MAX); +} + + +#endif // FFT_COMMON_HXX_ diff --git a/include/IRGridCommon.h b/include/IRGridCommon.h new file mode 100644 index 0000000..490eb37 --- /dev/null +++ b/include/IRGridCommon.h @@ -0,0 +1,1612 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : grid_common.hxx +// Author : Pavel A. Koshevoy +// Created : Wed Jan 10 09:31:00 MDT 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : code used to refine mesh transform control points + +#ifndef GRID_COMMON_HXX_ +#define GRID_COMMON_HXX_ + +// the includes: +#include "itkIRCommon.h" +#include "IRFFTCommon.h" +#include "IRGridTransform.h" +#include "itkGridTransform.h" +#include "itkRegularStepGradientDescentOptimizer2.h" +#include "IROptimizerObserver.h" +#include "IRLog.h" + +// system includes: +#include +#include +#include +#include +#include + +// ITK includes: +#include +#include +#include +#include +#include + +//---------------------------------------------------------------- +// DEBUG_REFINE_ONE_POINT +// +// #define DEBUG_REFINE_ONE_POINT + +//---------------------------------------------------------------- +// DEBUG_ESTIMATE_DISPLACEMENT +// +// #define DEBUG_ESTIMATE_DISPLACEMENT + +#if defined(DEBUG_REFINE_ONE_POINT) || defined(DEBUG_ESTIMATE_DISPLACEMENT) +static int COUNTER = 0; +#endif + + +//---------------------------------------------------------------- +// optimizer_t +// +typedef itk::RegularStepGradientDescentOptimizer2 optimizer_t; + + +//---------------------------------------------------------------- +// setup_grid_transform +// +extern bool +setup_grid_transform(the_grid_transform_t & transform, + unsigned int rows, + unsigned int cols, + const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const mask_t * tile_mask, + base_transform_t::ConstPointer mosaic_to_tile, + unsigned int max_iterations = 100, + double min_step_scale = 1e-12, + double min_error_sqrd = 1e-16, + unsigned int pick_up_pace_steps = 5); + +//---------------------------------------------------------------- +// setup_mesh_transform +// +extern bool +setup_mesh_transform(the_mesh_transform_t & transform, + unsigned int rows, + unsigned int cols, + const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const mask_t * tile_mask, + base_transform_t::ConstPointer mosaic_to_tile, + unsigned int max_iterations, + double min_step_scale, + double min_error_sqrd, + unsigned int pick_up_pace_steps); + +//---------------------------------------------------------------- +// estimate_displacement +// +template +void +estimate_displacement(the_log_t & log, + + // current best metric: + double & best_metric, + translate_transform_t::Pointer & best_transform, + + // origin offset of the small neighborhood of image A + // within the large image A: + const vec2d_t & offset, + + const TImage * a, // large image A + const TImage * b, // small image B + + // this could be the size of B (ususally when matching + // small neighborhoods in A and B), or the + // maximum dimesions of A and B (plain image matching): + const typename TImage::SizeType & period_sz, + + const local_max_t & lm, + const double overlap_min = 0.05, + const double overlap_max = 1.0, + const mask_t * mask_a = NULL, + const mask_t * mask_b = NULL, + + const unsigned int num_perms = 4) +{ + vec2d_t best_offset = best_transform->GetOffset(); + + const unsigned int & w = period_sz[0]; + const unsigned int & h = period_sz[1]; + + typedef typename TImage::SpacingType spacing_t; + spacing_t spacing = b->GetSpacing(); + double sx = spacing[0]; + double sy = spacing[1]; + + // evaluate 4 permutation of the maxima: + const double x = lm.x_; + const double y = lm.y_; + + const vec2d_t t[] = { + vec2d(sx * x, sy * y), vec2d(sx * (x - w), sy * y), vec2d(sx * x, sy * (y - h)), vec2d(sx * (x - w), sy * (y - h)) + }; + +#ifdef DEBUG_ESTIMATE_DISPLACEMENT + static const the_text_t fn_save("/tmp/estimate_displacement-"); + static int COUNTER2 = 0; + the_text_t suffix = the_text_t::number(COUNTER, 3, '0') + "-" + the_text_t::number(COUNTER2, 3, '0') + "-"; + COUNTER2++; + + // FIXME: + { + translate_transform_t::Pointer ti = translate_transform_t::New(); + ti->SetOffset(-offset); + save_rgb(fn_save + suffix + "0-init.png", a, b, ti, mask_a, mask_b); + } +#endif + + std::vector metricArray; + metricArray.reserve(num_perms); + + std::vector transformArray; + transformArray.reserve(num_perms); + + for (int i = 0; i < (int)num_perms; i++) + { + translate_transform_t::Pointer ti = translate_transform_t::New(); + ti->SetOffset(t[i] - offset); + + double area_ratio = overlap_ratio(a, b, ti); + log << i << ": " << setw(3) << int(area_ratio * 100.0) << "% of overlap, "; + + if (area_ratio < overlap_min || area_ratio > overlap_max) + { + log << "skipping..." << endl; + continue; + } + +#if 0 + typedef itk::LinearInterpolateImageFunction interpolator_t; + // typedef itk::MeanSquaresImageToImageMetric metric_t; + typedef itk::NormalizedCorrelationImageToImageMetric + metric_t; + double metric = eval_metric(ti, + a, + b, + mask_a, + mask_b); +#else + double metric = my_metric(area_ratio, a, b, ti, mask_a, mask_b, overlap_min, overlap_max); + metricArray.push_back(metric); +#endif + + transformArray.push_back(ti); + } + + for (int i = 0; i < (int)metricArray.size(); i++) + { + double metric = metricArray[i]; + translate_transform_t::Pointer ti = transformArray[i]; + + log << metric; + if (metric < best_metric) + { + best_offset = t[i]; + best_metric = metric; + log << " - better..." << endl; + +#ifdef DEBUG_ESTIMATE_DISPLACEMENT + // FIXME: + save_rgb(fn_save + suffix + the_text_t::number(i + 1) + "-perm.png", a, b, ti, mask_a, mask_b); +#endif + } + else + { + log << " - worse..." << endl; + } + } + + metricArray.clear(); + transformArray.clear(); + + best_transform->SetOffset(best_offset); +} + + +//---------------------------------------------------------------- +// match_one_pair +// +// match two images, find the best matching transform and +// return the corresponding match metric value +// +template +double +match_one_pair(the_log_t & log, + translate_transform_t::Pointer & best_transform, + + const TImage * fi, + const mask_t * fi_mask, + + const TImage * mi, + const mask_t * mi_mask, + + // this could be the size of B (ususally when matching + // small neighborhoods in A and B), or the + // maximum dimesions of A and B (plain image matching): + const typename TImage::SizeType & period_sz, + + const double & overlap_min, + const double & overlap_max, + + image_t::PointType offset_min, + image_t::PointType offset_max, + + const double & lp_filter_r, + const double & lp_filter_s, + + const unsigned int max_peaks, + const bool & consider_zero_displacement) +{ + static const vec2d_t offset = vec2d(0, 0); + best_transform = NULL; + + std::list max_list; + unsigned int total_peaks = + find_correlation(max_list, fi, mi, lp_filter_r, lp_filter_s, overlap_min, overlap_max); + + unsigned int num_peaks = reject_negligible_maxima(max_list, 2.0); + log << num_peaks << '/' << total_peaks << " eligible peaks, "; + + if (num_peaks > max_peaks) + { + log << "skipping..." << endl; + return std::numeric_limits::max(); + } + + double best_metric = std::numeric_limits::max(); + best_transform = translate_transform_t::New(); + best_transform->SetIdentity(); + + if (consider_zero_displacement) + { + // make sure to consider the no-displacement case: + estimate_displacement(log, + best_metric, + best_transform, + offset, + fi, + mi, + period_sz, + local_max_t(0, 0, 0, 0), + 0, + 1, + fi_mask, + mi_mask, + 1); // don't try permutations + } + + // choose the best peak: + for (std::list::iterator j = max_list.begin(); j != max_list.end(); ++j) + { + log << "evaluating permutations..." << endl; + const local_max_t & lm = *j; + + double prev_metric = best_metric; + vec2d_t prev_offset = best_transform->GetOffset(); + + estimate_displacement(log, + best_metric, + best_transform, + offset, + fi, + mi, + period_sz, + lm, // local maxima record + overlap_min, + overlap_max, + fi_mask, + mi_mask); + + // re-evaluate the overlap in the context of the small neighborhoods: + double overlap = overlap_ratio(fi, mi, best_transform); + if (overlap < overlap_min || overlap > overlap_max) + { + best_metric = prev_metric; + best_transform->SetOffset(prev_offset); + } + } + + return best_metric; +} + + +//---------------------------------------------------------------- +// match_one_pair +// +// match two images, find the best matching transform and +// return the corresponding match metric value +// +template +double +match_one_pair(the_log_t & log, + translate_transform_t::Pointer & best_transform, + + const TImage * fi_large, + const mask_t * fi_large_mask, + + const TImage * fi, + const mask_t * /* fi_mask */, + + const TImage * mi, + const mask_t * mi_mask, + + // origin offset of the small fixed image neighborhood + // within the full image: + const vec2d_t & offset, + + const double & overlap_min, + const double & overlap_max, + + // image_t::PointType offset_min, + // image_t::PointType offset_max, + + const double & lp_filter_r, + const double & lp_filter_s, + + const unsigned int max_peaks, + const bool & consider_zero_displacement) +{ + best_transform = NULL; + + std::list max_list; + unsigned int total_peaks = + find_correlation(max_list, fi, mi, lp_filter_r, lp_filter_s, overlap_min, overlap_max); + + unsigned int num_peaks = reject_negligible_maxima(max_list, 2.0); + log << num_peaks << '/' << total_peaks << " eligible peaks, "; + + if (num_peaks > max_peaks) + { + log << "skipping..." << endl; + return std::numeric_limits::max(); + } + + double best_metric = std::numeric_limits::max(); + best_transform = translate_transform_t::New(); + best_transform->SetIdentity(); + + typename TImage::SizeType period_sz = mi->GetLargestPossibleRegion().GetSize(); + + if (consider_zero_displacement) + { + // make sure to consider the no-displacement case: + estimate_displacement(log, + best_metric, + best_transform, + offset, + fi_large, + mi, + period_sz, + local_max_t(0, 0, 0, 0), + 0, + 1, + fi_large_mask, + mi_mask, + 1); // don't try permutations + } + + // choose the best peak: + for (std::list::iterator j = max_list.begin(); j != max_list.end(); ++j) + { + log << "evaluating permutations..." << endl; + const local_max_t & lm = *j; + + double prev_metric = best_metric; + vec2d_t prev_offset = best_transform->GetOffset(); + + estimate_displacement(log, + best_metric, + best_transform, + offset, + fi_large, + mi, + period_sz, + lm, // local maxima record + overlap_min, + overlap_max, + fi_large_mask, + mi_mask); + + // re-evaluate the overlap in the context of the small neighborhoods: + double overlap = overlap_ratio(fi, mi, best_transform); + if (overlap < overlap_min || overlap > overlap_max) + { + best_metric = prev_metric; + best_transform->SetOffset(prev_offset); + } + } + + return best_metric; +} + + +//---------------------------------------------------------------- +// refine_one_point_fft +// +template +bool +refine_one_point_fft(the_log_t & log, + + vec2d_t & shift, + const pnt2d_t & origin, + + // the large images and their masks: + const TImage * tile_0, + const mask_t * mask_0, + + // the extracted small neighborhoods and their masks: + const TImage * img_0, + const mask_t * msk_0, + const TImage * img_1, + const mask_t * msk_1) +{ + // feed the two neighborhoods into the FFT translation estimator: + translate_transform_t::Pointer translate; + double metric = match_one_pair(log, + + // best translation transform: + translate, + + // fixed image in full: + tile_0, + mask_0, + + // fixed image small neighborhood: + img_0, + msk_0, + + // moving image small neighborhood: + img_1, + msk_1, + + // origin offset of the small fixed image + // neighborhood in the large image: + vec2d(origin[0], origin[1]), + + 0.25, // overlap min + 1.0, // overlap max + 0.5, // low pass filter r + 0.1, // low pass filter s + 10, // number of peaks + true); // consider the no-displacement case + +#ifdef DEBUG_REFINE_ONE_POINT + static const the_text_t fn_save("/tmp/refine_one_point_fft-"); + the_text_t suffix = the_text_t::number(COUNTER, 3, '0'); + + save_composite(fn_save + suffix + "-a.png", img_0, img_1, identity_transform_t::New(), true); +#endif // DEBUG_REFINE_ONE_POINT + +#if defined(DEBUG_REFINE_ONE_POINT) || defined(DEBUG_ESTIMATE_DISPLACEMENT) + COUNTER++; +#endif + + if (metric == std::numeric_limits::max()) + { + return false; + } + +#ifdef DEBUG_REFINE_ONE_POINT + save_composite(fn_save + suffix + "-b.png", img_0, img_1, translate.GetPointer(), true); +#endif // DEBUG_REFINE_ONE_POINT + + shift = -translate->GetOffset(); + return true; +} + + +//---------------------------------------------------------------- +// refine_one_point_helper +// +template +bool +refine_one_point_helper(const pnt2d_t & center, + const pnt2d_t & origin, + const double & min_overlap, + + // the large images and their masks: + const TImage * tile_0, + const mask_t * mask_0, + const TImage * tile_1, + const mask_t * mask_1, + + // transform that maps from the space of + // tile 0 to the space of tile 1: + const base_transform_t * t_01, + + // the extracted small neighborhoods and their masks: + TImage * img_0, + mask_t * msk_0, + TImage * img_1, + mask_t * msk_1, + + // neighborhood size and pixel spacing: + const typename TImage::SizeType & sz, + const typename TImage::SpacingType & sp) +{ + // setup the image interpolators: + typedef itk::LinearInterpolateImageFunction interpolator_t; + typename interpolator_t::Pointer interpolator[] = { interpolator_t::New(), interpolator_t::New() }; + + interpolator[0]->SetInputImage(tile_0); + interpolator[1]->SetInputImage(tile_1); + + if (!interpolator[0]->IsInsideBuffer(center)) + { + // don't bother extracting neigborhoods if the neighborhood center + // doesn't fall inside both images: + return false; + } + + // temporaries: + pnt2d_t mosaic_pt; + pnt2d_t tile_pt; + typename TImage::IndexType index; + + // keep count of pixel in the neighborhood: + unsigned int area[] = { 0, 0 }; + + // extract a neighborhood of the given point from both tiles: + for (unsigned int y = 0; y < sz[1]; y++) + { + mosaic_pt[1] = origin[1] + double(y) * sp[1]; + index[1] = y; + for (unsigned int x = 0; x < sz[0]; x++) + { + mosaic_pt[0] = origin[0] + double(x) * sp[0]; + index[0] = x; + + // fixed image: + tile_pt = mosaic_pt; + if (interpolator[0]->IsInsideBuffer(tile_pt) && pixel_in_mask(mask_0, tile_pt)) + { + double p = interpolator[0]->Evaluate(tile_pt); + img_0->SetPixel(index, (unsigned char)(std::min(255.0, p))); + msk_0->SetPixel(index, 1); + area[0]++; + } + else + { + img_0->SetPixel(index, 0); + msk_0->SetPixel(index, 0); + } + + // moving image: + tile_pt = t_01->TransformPoint(mosaic_pt); + if (interpolator[1]->IsInsideBuffer(tile_pt) && pixel_in_mask(mask_1, tile_pt)) + { + double p = interpolator[1]->Evaluate(tile_pt); + img_1->SetPixel(index, (unsigned char)(std::min(255.0, p))); + msk_1->SetPixel(index, 1); + area[1]++; + } + else + { + img_1->SetPixel(index, 0); + msk_1->SetPixel(index, 0); + } + } + } + + // skip points which don't have enough neighborhood information: + double max_area = double(sz[0] * sz[1]); + double a[] = { double(area[0]) / max_area, double(area[1]) / max_area }; + + if (a[0] < min_overlap || a[1] < min_overlap) + { + return false; + } + + return true; +} + +//---------------------------------------------------------------- +// refine_one_point_fft +// +template +bool +refine_one_point_fft(the_log_t & log, + vec2d_t & shift, + const pnt2d_t & center, + const pnt2d_t & origin, + const double & min_overlap, + + // the large images and their masks: + const TImage * tile_0, + const mask_t * mask_0, + const TImage * tile_1, + const mask_t * mask_1, + + // transform that maps from the space of + // tile 0 to the space of tile 1: + const base_transform_t * t_01, + + // the extracted small neighborhoods and their masks: + TImage * img_0, + mask_t * msk_0, + TImage * img_1, + mask_t * msk_1, + + // neighborhood size and pixel spacing: + const typename TImage::SizeType & sz, + const typename TImage::SpacingType & sp) +{ + if (!refine_one_point_helper( + center, origin, min_overlap, tile_0, mask_0, tile_1, mask_1, t_01, img_0, msk_0, img_1, msk_1, sz, sp)) + { + return false; + } + + return refine_one_point_fft(log, shift, origin, tile_0, mask_0, img_0, msk_0, img_1, msk_1); +} + +//---------------------------------------------------------------- +// refine_one_point_fft +// +template +bool +refine_one_point_fft(the_log_t & log, + + vec2d_t & shift, + const pnt2d_t & center, + const double & min_overlap, + + const TImage * tile_0, + const mask_t * mask_0, + const TImage * tile_1, + const mask_t * mask_1, + + // transform that maps from the space of + // tile 0 to the space of tile 1: + const base_transform_t * t_01, + + // neighborhood size: + const unsigned int & neighborhood) +{ + typename TImage::SpacingType sp = tile_1->GetSpacing(); + typename TImage::SizeType sz; + sz[0] = neighborhood; + sz[1] = neighborhood; + + pnt2d_t origin(center); + origin[0] -= 0.5 * double(sz[0]) * sp[0]; + origin[1] -= 0.5 * double(sz[1]) * sp[1]; + + typename TImage::Pointer img[] = { make_image(sz), make_image(sz) }; + + mask_t::Pointer msk[] = { make_image(sz), make_image(sz) }; + + img[0]->SetSpacing(sp); + img[1]->SetSpacing(sp); + msk[0]->SetSpacing(sp); + msk[1]->SetSpacing(sp); + + return refine_one_point_fft(log, + shift, + center, + origin, + min_overlap, + tile_0, + mask_0, + tile_1, + mask_1, + t_01, + img[0], + msk[0], + img[1], + msk[1], + sz, + sp); +} + +//---------------------------------------------------------------- +// refine_one_point_fft +// +template +bool +refine_one_point_fft(the_log_t & log, + vec2d_t & shift, + + // the extracted small neighborhoods and their masks: + const TImage * img_0, + const mask_t * msk_0, + const TImage * img_1, + const mask_t * msk_1) +{ + typename TImage::SizeType period_sz = img_1->GetLargestPossibleRegion().GetSize(); + + // feed the two neighborhoods into the FFT translation estimator: + translate_transform_t::Pointer translate; + double metric = match_one_pair(log, + + // best translation transform: + translate, + + // fixed image small neighborhood: + img_0, + msk_0, + + // moving image small neighborhood: + img_1, + msk_1, + + period_sz, + + 0.25, // overlap min + 1.0, // overlap max + 0.5, // low pass filter r + 0.1, // low pass filter s + 10, // number of peaks + true); // consider the no-displacement case + +#ifdef DEBUG_REFINE_ONE_POINT + static const the_text_t fn_save("/tmp/refine_one_point_fft-"); + the_text_t suffix = the_text_t::number(COUNTER, 3, '0'); + + save_composite(fn_save + suffix + "-a.png", img_0, img_1, identity_transform_t::New(), true); +#endif // DEBUG_REFINE_ONE_POINT + +#if defined(DEBUG_REFINE_ONE_POINT) || defined(DEBUG_ESTIMATE_DISPLACEMENT) + COUNTER++; +#endif + + if (metric == std::numeric_limits::max()) + { + return false; + } + +#ifdef DEBUG_REFINE_ONE_POINT + save_composite(fn_save + suffix + "-b.png", img_0, img_1, translate.GetPointer(), true); +#endif // DEBUG_REFINE_ONE_POINT + + shift = -translate->GetOffset(); + return true; +} + + +//---------------------------------------------------------------- +// extract +// +// Return the number of pixels in the region in the mask: +// +template +unsigned int +extract(const typename TImage::PointType & origin, + + const TImage * tile, + const TMask * mask, + + TImage * tile_region, + TMask * mask_region, + + const typename TImage::PixelType & bg = itk::NumericTraits::Zero) +{ + tile_region->SetOrigin(origin); + mask_region->SetOrigin(origin); + + // setup the image interpolator: + typedef itk::LinearInterpolateImageFunction interpolator_t; + typename interpolator_t::Pointer interpolator = interpolator_t::New(); + interpolator->SetInputImage(tile); + + // temporary: + typename TImage::PointType tile_pt; + typename TImage::PixelType tile_val; + typename TMask::IndexType mask_ix; + typename TMask::PixelType mask_val; + + static const typename TMask::PixelType mask_min = itk::NumericTraits::Zero; + + static const typename TMask::PixelType mask_max = itk::NumericTraits::One; + + unsigned int num_pixels = 0; + + typedef itk::ImageRegionIteratorWithIndex itex_t; + itex_t itex(tile_region, tile_region->GetLargestPossibleRegion()); + for (itex.GoToBegin(); !itex.IsAtEnd(); ++itex) + { + tile_val = bg; + mask_val = mask_min; + + tile_region->TransformIndexToPhysicalPoint(itex.GetIndex(), tile_pt); + if (interpolator->IsInsideBuffer(tile_pt)) + { + mask_val = mask_max; + if (mask != NULL) + { + mask->TransformPhysicalPointToIndex(tile_pt, mask_ix); + mask_val = mask->GetPixel(mask_ix); + } + + if (mask_val != mask_min) + { + tile_val = (typename TImage::PixelType)(interpolator->Evaluate(tile_pt)); + num_pixels++; + } + } + + itex.Set(tile_val); + mask_region->SetPixel(itex.GetIndex(), mask_val); + } + + return num_pixels; +} + + +//---------------------------------------------------------------- +// refine_one_point_helper +// +template +bool +refine_one_point_helper( // the large images and their masks: + const TImage * tile_0, + const mask_t * mask_0, + const TImage * tile_1, + const mask_t * mask_1, + + + // mosaic space neighborhood center: + const pnt2d_t & center, + pnt2d_t & origin, + + // minimum acceptable neighborhood overlap ratio: + const double & min_overlap, + + // the extracted small neighborhoods and their masks: + TImage * img_0, + mask_t * msk_0, + TImage * img_1, + mask_t * msk_1) +{ + const typename TImage::SizeType & sz = img_0->GetLargestPossibleRegion().GetSize(); + const typename TImage::SpacingType & sp = img_0->GetSpacing(); + + origin[0] = center[0] - (double(sz[0]) * sp[0]) / 2; + origin[1] = center[1] - (double(sz[1]) * sp[1]) / 2; + + pnt2d_t corner[] = { origin + vec2d(sz[0] * sp[0], 0), + origin + vec2d(sz[0] * sp[0], sz[1] * sp[1]), + origin + vec2d(0, sz[1] * sp[1]) }; + + // make sure both tiles overlap the region: + typename TImage::IndexType tmp_ix; + if (!(tile_0->TransformPhysicalPointToIndex(origin, tmp_ix) || + tile_0->TransformPhysicalPointToIndex(corner[0], tmp_ix) || + tile_0->TransformPhysicalPointToIndex(corner[1], tmp_ix) || + tile_0->TransformPhysicalPointToIndex(corner[2], tmp_ix)) || + !(tile_1->TransformPhysicalPointToIndex(origin, tmp_ix) || + tile_1->TransformPhysicalPointToIndex(corner[0], tmp_ix) || + tile_1->TransformPhysicalPointToIndex(corner[1], tmp_ix) || + tile_1->TransformPhysicalPointToIndex(corner[2], tmp_ix))) + { + return false; + } + + // keep count of pixel in the neighborhood: + unsigned int area[] = { 0, 0 }; + + area[0] = extract(origin, tile_0, mask_0, img_0, msk_0); + + area[1] = extract(origin, tile_1, mask_1, img_1, msk_1); + + // skip points which don't have enough neighborhood information: + double max_area = double(sz[0] * sz[1]); + double a[] = { double(area[0]) / max_area, double(area[1]) / max_area }; + + if (a[0] < min_overlap || a[1] < min_overlap) + { + return false; + } + + return true; +} + +//---------------------------------------------------------------- +// tiles_intersect +// +template +bool +tiles_intersect( // the large images and their masks: + const TImage * tile_0, + const TImage * tile_1, + + // transforms from mosaic space to tile space: + const base_transform_t * forward_0, + const base_transform_t * forward_1, + + // mosaic space neighborhood center: + const pnt2d_t & origin, + + // neighborhood size and pixel spacing: + const typename TImage::SizeType & sz, + const typename TImage::SpacingType & sp) +{ + pnt2d_t corner[] = { origin, + origin + vec2d(sz[0] * sp[0], 0), + origin + vec2d(sz[0] * sp[0], sz[1] * sp[1]), + origin + vec2d(0, sz[1] * sp[1]) }; + + // temporary: + typename TImage::IndexType tmp_ix; + pnt2d_t tmp_pt; + + // make sure both tiles overlap the region: + bool in[] = { false, false }; + + for (unsigned int i = 0; i < 4 && (!in[0] || !in[1]); i++) + { + tmp_pt = forward_0->TransformPoint(corner[i]); + in[0] = in[0] || tile_0->TransformPhysicalPointToIndex(tmp_pt, tmp_ix); + + tmp_pt = forward_1->TransformPoint(corner[i]); + in[1] = in[1] || tile_1->TransformPhysicalPointToIndex(tmp_pt, tmp_ix); + } + + return in[0] && in[1]; +} + +//---------------------------------------------------------------- +// refine_one_point_helper +// +template +bool +refine_one_point_helper( // the large images and their masks: + const TImage * tile_0, + const mask_t * mask_0, + const TImage * tile_1, + const mask_t * mask_1, + + // transforms from mosaic space to tile space: + const base_transform_t * forward_0, + const base_transform_t * forward_1, + + // mosaic space neighborhood center: + const pnt2d_t & center, + + // minimum acceptable neighborhood overlap ratio: + const double & min_overlap, + + // neighborhood size and pixel spacing: + const typename TImage::SizeType & sz, + const typename TImage::SpacingType & sp, + + // the extracted small neighborhoods and their masks: + TImage * img_0_large, + mask_t * msk_0_large, + TImage * img_0, + mask_t * msk_0, + TImage * img_1, + mask_t * msk_1) +{ + pnt2d_t origin = center; + origin[0] -= (double(sz[0]) * sp[0]) / 2; + origin[1] -= (double(sz[1]) * sp[1]) / 2; + + if (!tiles_intersect(tile_0, tile_1, forward_0, forward_1, origin, sz, sp)) + { + return false; + } + + // setup the image interpolators: + typedef itk::LinearInterpolateImageFunction interpolator_t; + typename interpolator_t::Pointer interpolator[] = { interpolator_t::New(), interpolator_t::New() }; + + interpolator[0]->SetInputImage(tile_0); + interpolator[1]->SetInputImage(tile_1); + + // temporaries: + pnt2d_t mosaic_pt; + pnt2d_t tile_pt; + typename TImage::IndexType index; + + // keep count of pixel in the neighborhood: + unsigned int area[] = { 0, 0 }; + + // extract a neighborhood of the given point from both tiles: + for (unsigned int y = 0; y < sz[1]; y++) + { + mosaic_pt[1] = origin[1] + double(y) * sp[1]; + index[1] = y; + for (unsigned int x = 0; x < sz[0]; x++) + { + mosaic_pt[0] = origin[0] + double(x) * sp[0]; + index[0] = x; + + // fixed image: + tile_pt = forward_0->TransformPoint(mosaic_pt); + if (interpolator[0]->IsInsideBuffer(tile_pt) && pixel_in_mask(mask_0, tile_pt)) + { + double p = interpolator[0]->Evaluate(tile_pt); + img_0->SetPixel(index, (unsigned char)(std::min(255.0, p))); + msk_0->SetPixel(index, 1); + area[0]++; + } + else + { + img_0->SetPixel(index, 0); + msk_0->SetPixel(index, 0); + } + + // moving image: + tile_pt = forward_1->TransformPoint(mosaic_pt); + if (interpolator[1]->IsInsideBuffer(tile_pt) && pixel_in_mask(mask_1, tile_pt)) + { + double p = interpolator[1]->Evaluate(tile_pt); + img_1->SetPixel(index, (unsigned char)(std::min(255.0, p))); + msk_1->SetPixel(index, 1); + area[1]++; + } + else + { + img_1->SetPixel(index, 0); + msk_1->SetPixel(index, 0); + } + } + } + + // skip points which don't have enough neighborhood information: + double max_area = double(sz[0] * sz[1]); + double a[] = { double(area[0]) / max_area, double(area[1]) / max_area }; + + if (a[0] < min_overlap || a[1] < min_overlap) + { + return false; + } + + if (img_0_large != NULL) + { + // extract the larger neighborhood from the fixed tile: + pnt2d_t origin_large(center); + origin_large[0] -= sz[0] * sp[0]; + origin_large[1] -= sz[1] * sp[1]; + + for (unsigned int y = 0; y < sz[1] * 2; y++) + { + mosaic_pt[1] = origin_large[1] + double(y) * sp[1]; + index[1] = y; + for (unsigned int x = 0; x < sz[0] * 2; x++) + { + mosaic_pt[0] = origin_large[0] + double(x) * sp[0]; + index[0] = x; + + // fixed image: + tile_pt = forward_0->TransformPoint(mosaic_pt); + if (interpolator[0]->IsInsideBuffer(tile_pt) && pixel_in_mask(mask_0, tile_pt)) + { + double p = interpolator[0]->Evaluate(tile_pt); + img_0_large->SetPixel(index, (unsigned char)(std::min(255.0, p))); + msk_0_large->SetPixel(index, 1); + } + else + { + img_0_large->SetPixel(index, 0); + msk_0_large->SetPixel(index, 0); + } + } + } + } + + return true; +} + + +//---------------------------------------------------------------- +// refine_one_point_fft +// +template +bool +refine_one_point_fft(the_log_t & log, + vec2d_t & shift, + + // fixed tile, in mosaic space: + const TImage * tile_0, + const mask_t * mask_0, + + // moving tile, in mosaic space: + const TImage * tile_1, + const mask_t * mask_1, + + // mosaic space neighborhood center: + const pnt2d_t & center, + + // minimum acceptable neighborhood overlap ratio: + const double & min_overlap, + + // the extracted small neighborhoods and their masks: + TImage * img_0, + mask_t * msk_0, + TImage * img_1, + mask_t * msk_1) +{ + pnt2d_t origin; + if (!refine_one_point_helper( + tile_0, mask_0, tile_1, mask_1, center, origin, min_overlap, img_0, msk_0, img_1, msk_1)) + { + return false; + } + + return refine_one_point_fft(log, + shift, + pnt2d(0, 0), // origin, + tile_0, + mask_0, + img_0, + msk_0, + img_1, + msk_1); +} + + +//---------------------------------------------------------------- +// refine_one_point_fft +// +template +bool +refine_one_point_fft(the_log_t & log, + vec2d_t & shift, + + // the large images and their masks: + const TImage * tile_0, + const mask_t * mask_0, + const TImage * tile_1, + const mask_t * mask_1, + + // transforms from mosaic space to tile space: + const base_transform_t * forward_0, + const base_transform_t * forward_1, + + // mosaic space neighborhood center: + const pnt2d_t & center, + + // minimum acceptable neighborhood overlap ratio: + const double & min_overlap, + + // neighborhood size and pixel spacing: + const typename TImage::SizeType & sz, + const typename TImage::SpacingType & sp, + + // the extracted larger neighborhood: + TImage * img_0_large, + mask_t * msk_0_large, + + // the extracted small neighborhoods and their masks: + TImage * img_0, + mask_t * msk_0, + TImage * img_1, + mask_t * msk_1) +{ + if (!refine_one_point_helper(tile_0, + mask_0, + tile_1, + mask_1, + forward_0, + forward_1, + center, + min_overlap, + sz, + sp, + img_0_large, + msk_0_large, + img_0, + msk_0, + img_1, + msk_1)) + { + return false; + } + + pnt2d_t origin; + origin[0] = sz[0] * sp[0] / 2; + origin[1] = sz[1] * sp[1] / 2; + return refine_one_point_fft(log, shift, origin, img_0_large, msk_0_large, img_0, msk_0, img_1, msk_1); +} + + +//---------------------------------------------------------------- +// refine_one_point_fft +// +template +bool +refine_one_point_fft(the_log_t & log, + vec2d_t & shift, + + // the large images and their masks: + const TImage * tile_0, + const mask_t * mask_0, + const TImage * tile_1, + const mask_t * mask_1, + + // transforms from mosaic space to tile space: + const base_transform_t * forward_0, + const base_transform_t * forward_1, + + // mosaic space neighborhood center: + const pnt2d_t & center, + + // minimum acceptable neighborhood overlap ratio: + const double & min_overlap, + + // neighborhood size and pixel spacing: + const typename TImage::SizeType & sz, + const typename TImage::SpacingType & sp, + + // the extracted small neighborhoods and their masks: + TImage * img_0, + mask_t * msk_0, + TImage * img_1, + mask_t * msk_1) +{ + pnt2d_t origin(center); + origin[0] -= sz[0] * sp[0] / 2; + origin[1] -= sz[1] * sp[1] / 2; + + if (!refine_one_point_helper(tile_0, + mask_0, + tile_1, + mask_1, + forward_0, + forward_1, + origin, + min_overlap, + sz, + sp, + NULL, + NULL, + img_0, + msk_0, + img_1, + msk_1)) + { + return false; + } + + return refine_one_point_fft(log, shift, img_0, msk_0, img_1, msk_1); +} + + +//---------------------------------------------------------------- +// refine_one_pair +// +template +double +refine_one_pair(the_log_t & log, + + typename TTransform::Pointer & t01, + + const TImage * i0, + const mask_t * m0, + + const TImage * i1, + const mask_t * m1, + + const unsigned int levels, + const unsigned int iterations, + const double & min_step, + const double & max_step) +{ + // typedef itk::MultiResolutionImageRegistrationMethod + typedef itk::ImageRegistrationMethod registration_t; + + // setup the registration object: + typename registration_t::Pointer registration = registration_t::New(); + + // setup the image interpolator: + typedef itk::LinearInterpolateImageFunction interpolator_t; + registration->SetInterpolator(interpolator_t::New()); + + // setup the optimizer: + optimizer_t::Pointer optimizer = optimizer_t::New(); + registration->SetOptimizer(optimizer); + optimizer_observer_t::Pointer observer = optimizer_observer_t::New(); + observer->log_ = &log; + optimizer->SetLog(&log); + optimizer->AddObserver(itk::IterationEvent(), observer); + optimizer->SetMinimize(true); + optimizer->SetNumberOfIterations(iterations); + optimizer->SetMinimumStepLength(min_step); + optimizer->SetMaximumStepLength(max_step); + optimizer->SetGradientMagnitudeTolerance(1e-6); + optimizer->SetRelaxationFactor(5e-1); + optimizer->SetPickUpPaceSteps(5); + optimizer->SetBackTracking(true); + + // FIXME: this is probably unnecessary: + typedef optimizer_t::ScalesType optimizer_scales_t; + optimizer_scales_t parameter_scales(t01->GetNumberOfParameters()); + parameter_scales.Fill(1.0); + + try + { + optimizer->SetScales(parameter_scales); + } + catch (itk::ExceptionObject &) + {} + + // setup the image-to-image metric: + typedef itk::NormalizedCorrelationImageToImageMetric metric_t; + + typename metric_t::Pointer metric = metric_t::New(); + registration->SetMetric(metric); + + registration->SetTransform(t01); + registration->SetInitialTransformParameters(t01->GetParameters()); + + // setup the masks: + typedef itk::ImageMaskSpatialObject<2> mask_so_t; + mask_t::ConstPointer fi_mask = m0; + if (m0 != NULL) + { + mask_so_t::Pointer fi_mask_so = mask_so_t::New(); + fi_mask_so->SetImage(fi_mask); + metric->SetFixedImageMask(fi_mask_so); + } + + mask_t::ConstPointer mi_mask = m1; + if (m1 != NULL) + { + mask_so_t::Pointer mi_mask_so = mask_so_t::New(); + mi_mask_so->SetImage(mi_mask); + metric->SetMovingImageMask(mi_mask_so); + } + + // setup the fixed and moving image: + typename TImage::ConstPointer fi = i0; + typename TImage::ConstPointer mi = i1; + + registration->SetFixedImageRegion(fi->GetLargestPossibleRegion()); + registration->SetFixedImage(fi); + registration->SetMovingImage(mi); + + // setup the multi-resolution image pyramids: + // registration->SetNumberOfLevels(levels); + + // evaluate the metric before the registration: + double metric_before = eval_metric(t01, fi, mi, fi_mask, mi_mask); + assert(metric_before != std::numeric_limits::max()); + + typename TTransform::ParametersType params_before = t01->GetParameters(); + + // perform the registration: + try + { + log << endl; + registration->StartRegistration(); + } + catch (itk::ExceptionObject & exception) + { + log << "image registration threw an exception: " << endl << exception.what() << endl; + } + t01->SetParameters(optimizer->GetBestParams()); + + // evaluate the metric after the registration: + double metric_after = eval_metric(t01, fi, mi, fi_mask, mi_mask); + + typename TTransform::ParametersType params_after = t01->GetParameters(); + + log << "BEFORE: " << metric_before << endl << "AFTER: " << metric_after << endl; + + if (metric_before <= metric_after || metric_after != metric_after) + { + log << "NOTE: minimization failed, ignoring registration results..." << endl; + t01->SetParameters(params_before); + return metric_before; + } + + return metric_after; +} + + +//---------------------------------------------------------------- +// refine_one_point_nofft +// +template +bool +refine_one_point_nofft(the_log_t & log, + + vec2d_t & shift, + const pnt2d_t & origin, + + // the large images and their masks: + const TImage * tile_0, + const mask_t * mask_0, + + // the extracted small neighborhood with masks: + const TImage * img_1, + const mask_t * msk_1) +{ + // feed the two neighborhoods into the optimizer translation estimator: + translate_transform_t::Pointer translate = translate_transform_t::New(); + vec2d_t offset = vec2d(origin[0], origin[1]); + translate->SetOffset(-offset); + + double metric = refine_one_pair(translate, + tile_0, + mask_0, + img_1, + msk_1, + 3, // number of pyramid levels + 100, // number of iterations + 1e-8, // min step + 1e+4); // max step + + shift = -(translate->GetOffset() + offset); + log << "shift: " << shift << endl; + +#ifdef DEBUG_REFINE_ONE_POINT + static const the_text_t fn_save("/tmp/refine_one_point_nofft-"); + the_text_t suffix = the_text_t::number(COUNTER, 3, '0'); + { + translate_transform_t::Pointer t = translate_transform_t::New(); + t->SetOffset(-offset); + save_composite(fn_save + suffix + "-a.png", tile_0, img_1, t.GetPointer(), true); + } +#endif // DEBUG_REFINE_ONE_POINT + +#if defined(DEBUG_REFINE_ONE_POINT) || defined(DEBUG_ESTIMATE_DISPLACEMENT) + COUNTER++; +#endif + + if (metric == std::numeric_limits::max()) + { + return false; + } + +#ifdef DEBUG_REFINE_ONE_POINT + save_composite(fn_save + suffix + "-b.png", tile_0, img_1, translate.GetPointer(), true); +#endif // DEBUG_REFINE_ONE_POINT + + return true; +} + + +//---------------------------------------------------------------- +// refine_one_point_nofft +// +template +bool +refine_one_point_nofft(vec2d_t & shift, + const pnt2d_t & center, + const pnt2d_t & origin, + const double & min_overlap, + + // the large images and their masks: + const TImage * tile_0, + const mask_t * mask_0, + const TImage * tile_1, + const mask_t * mask_1, + + // transform that maps from the space of + // tile 0 to the space of tile 1: + const base_transform_t * t_01, + + // the extracted small neighborhoods and their masks: + TImage * img_0, + mask_t * msk_0, + TImage * img_1, + mask_t * msk_1, + + // neighborhood size and pixel spacing: + const typename TImage::SizeType & sz, + const typename TImage::SpacingType & sp) +{ + if (!refine_one_point_helper( + center, origin, min_overlap, tile_0, mask_0, tile_1, mask_1, t_01, img_0, msk_0, img_1, msk_1, sz, sp)) + { + return false; + } + + // feed the two neighborhoods into the optimizer translation estimator: + translate_transform_t::Pointer translate = translate_transform_t::New(); + translate->SetIdentity(); + + double metric = refine_one_pair(translate, + img_0, + msk_0, + img_1, + msk_1, + 3, // number of pyramid levels + 100, // number of iterations + 1e-8, // min step + 1e+4); // max step + +#ifdef DEBUG_REFINE_ONE_POINT + static const the_text_t fn_save("/tmp/refine_one_point_nofft-"); + the_text_t suffix = the_text_t::number(COUNTER, 3, '0'); + save_composite( + fn_save + suffix + "-a.png", img_0.GetPointer(), img_1.GetPointer(), identity_transform_t::New(), true); +#endif // DEBUG_REFINE_ONE_POINT + +#if defined(DEBUG_REFINE_ONE_POINT) || defined(DEBUG_ESTIMATE_DISPLACEMENT) + COUNTER++; +#endif + + if (metric == std::numeric_limits::max()) + { + return false; + } + +#ifdef DEBUG_REFINE_ONE_POINT + save_composite(fn_save + suffix + "-b.png", img_0.GetPointer(), img_1.GetPointer(), translate.GetPointer(), true); +#endif // DEBUG_REFINE_ONE_POINT + + shift = -translate->GetOffset(); + return true; +} + +//---------------------------------------------------------------- +// refine_one_point_nofft +// +template +bool +refine_one_point_nofft(vec2d_t & shift, + const pnt2d_t & center, + const double & min_overlap, + const TImage * tile_0, + const mask_t * mask_0, + const TImage * tile_1, + const mask_t * mask_1, + + // transform that maps from the space of + // tile 0 to the space of tile 1: + const base_transform_t * t_01, + + // neighborhood size: + const unsigned int & neighborhood) +{ + typename TImage::SpacingType sp = tile_1->GetSpacing(); + typename TImage::SizeType sz; + sz[0] = neighborhood; + sz[1] = neighborhood; + + pnt2d_t origin(center); + origin[0] -= 0.5 * double(sz[0]) * sp[0]; + origin[1] -= 0.5 * double(sz[1]) * sp[1]; + + typename TImage::Pointer img[] = { make_image(sz), make_image(sz) }; + + mask_t::Pointer msk[] = { make_image(sz), make_image(sz) }; + + img[0]->SetSpacing(sp); + img[1]->SetSpacing(sp); + msk[0]->SetSpacing(sp); + msk[1]->SetSpacing(sp); + + return refine_one_point_nofft( + shift, center, origin, min_overlap, tile_0, mask_0, tile_1, mask_1, t_01, img[0], msk[0], img[1], msk[1], sz, sp); +} + + +#endif // GRID_COMMON_HXX_ diff --git a/include/IRGridTransform.h b/include/IRGridTransform.h new file mode 100644 index 0000000..b4b0099 --- /dev/null +++ b/include/IRGridTransform.h @@ -0,0 +1,272 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: nil -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_grid_transform.hxx +// Author : Pavel A. Koshevoy +// Created : Thu Nov 30 10:45:47 MST 2006 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A discontinuous transform -- a uniform grid of vertices is +// mapped to an image. At each vertex, in addition to image +// space coordinates, a second set of coordinates is stored. +// This is similar to texture mapped OpenGL triangle meshes, +// where the texture coordinates correspond to the image space +// vertex coordinates. + +#ifndef THE_GRID_TRANSFORM_HXX_ +#define THE_GRID_TRANSFORM_HXX_ + +// local includes: +#include "itkIRCommon.h" + +// system includes: +#include +#include + + +//---------------------------------------------------------------- +// vertex_t +// +class vertex_t +{ +public: + // normalized tile space coordinates, typically [0, 1] x [0, 1]: + pnt2d_t uv_; + + // physical space coordinates: + pnt2d_t xy_; +}; + +//---------------------------------------------------------------- +// triangle_t +// +class triangle_t +{ +public: + triangle_t(); + + // check whether a given xy-point falls within this triangle, + // return corresponding uv-point (not triangle barycentric coordinates): + bool + xy_intersect(const vertex_t * v_arr, const pnt2d_t & xy, pnt2d_t & uv) const; + + // check whether a given uv-point falls within this triangle, + // return corresponding xy-point (not triangle barycentric coordinates): + bool + uv_intersect(const vertex_t * v_arr, const pnt2d_t & uv, pnt2d_t & xy) const; + + // triangle vertex indices, counterclockwise winding: + unsigned int vertex_[3]; + + // precomputed fast barycentric coordinate calculation coefficients: + + // for intersection calculation in xy-space: + double xy_pwb[3]; + double xy_pwc[3]; + + // for intersection calculation in uv-space: + double uv_pwb[3]; + double uv_pwc[3]; +}; + +//---------------------------------------------------------------- +// the_acceleration_grid_t +// +// The bounding grid triangle/point intersection acceleration +// structure used to speed up grid transform and mesh transform: +// +class the_acceleration_grid_t +{ +public: + the_acceleration_grid_t(); + + // find the grid cell containing a given xy-point: + unsigned int + xy_cell(const pnt2d_t & xy) const; + + // find the triangle containing a given xy-point, + // and calculate corresponding uv-point: + unsigned int + xy_triangle(const pnt2d_t & xy, pnt2d_t & uv) const; + + // find the grid cell containing a given uv-point: + unsigned int + uv_cell(const pnt2d_t & uv) const; + + // find the triangle containing a given uv-point, + // and calculate corresponding xy-point: + unsigned int + uv_triangle(const pnt2d_t & uv, pnt2d_t & xy) const; + + // update the vertex xy coordinates and rebuild the grid + void + update(const vec2d_t * xy_shift); + void + shift(const vec2d_t & xy_shift); + + // resize the grid: + void + resize(unsigned int rows, unsigned int cols); + + // rebuild the acceleration grid: + void + rebuild(); + +private: + // helper used to rebuild the grid: + void + update_grid(unsigned int t_idx); + +public: + // the acceleration structure: + std::vector> xy_; + std::vector> uv_; + unsigned int rows_; + unsigned int cols_; + + // the grid bounding box (in xy-space): + pnt2d_t xy_min_; + vec2d_t xy_ext_; + + // the triangle mesh: + std::vector mesh_; + std::vector tri_; +}; + + +//---------------------------------------------------------------- +// the_base_triangle_transform_t +// +class the_base_triangle_transform_t +{ +public: + the_base_triangle_transform_t() {} + + // transform the point: + bool + transform(const pnt2d_t & xy, pnt2d_t & uv) const; + + // inverse transform the point: + bool + transform_inv(const pnt2d_t & uv, pnt2d_t & xy) const; + + // calculate the derivatives of the transforms with respect to + // transform parameters: + bool + jacobian(const pnt2d_t & xy, unsigned int * idx, double * jac) const; + +public: + // tile bounding box: + pnt2d_t tile_min_; + vec2d_t tile_ext_; + + // the acceleration grid (stores triangle vertices, and triangles): + the_acceleration_grid_t grid_; +}; + + +//---------------------------------------------------------------- +// the_grid_transform_t +// +class the_grid_transform_t : public the_base_triangle_transform_t +{ +public: + the_grid_transform_t(); + + // check to see whether the transform has already been setup: + bool + is_ready() const; + + // vertex accessors: + inline const vertex_t & + vertex(size_t row, size_t col) const + { + return grid_.mesh_[row * (cols_ + 1) + col]; + } + + inline vertex_t & + vertex(size_t row, size_t col) + { + return grid_.mesh_[row * (cols_ + 1) + col]; + } + + // inverse transform the point: + bool + transform_inv(const pnt2d_t & uv, pnt2d_t & xy) const; + + // setup the transform: + void + setup(unsigned int rows, + unsigned int cols, + const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const std::vector & xy); + +private: + // helper used to setup the triangle mesh: + void + setup_mesh(); + +public: + // number of rows and columns of quads in the mesh + // (each quad is made up of 2 triangles): + size_t rows_; + size_t cols_; +}; + + +//---------------------------------------------------------------- +// the_mesh_transform_t +// +class the_mesh_transform_t : public the_base_triangle_transform_t +{ +public: + // check to see whether the transform has already been setup: + bool + is_ready() const; + + // setup the transform: + bool + setup(const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const std::vector & uv, + const std::vector & xy, + unsigned int accel_grid_rows = 16, + unsigned int accel_grid_cols = 16); + + // insert a point into the mesh, and re-triangulate + // using Delaunay triangulation: + bool + insert_point(const pnt2d_t & uv, const pnt2d_t & xy, const bool delay_setup = false); + + // insert a point into the mesh (xy-point is extrapolated), + // and re-triangulate using Delaunay triangulation: + bool + insert_point(const pnt2d_t & uv); + +private: + // helper used to setup the triangle mesh: + bool + setup_mesh(); +}; + + +#endif // THE_GRID_TRANSFORM_HXX_ diff --git a/include/IRMosaicRefinementCommon.h b/include/IRMosaicRefinementCommon.h new file mode 100644 index 0000000..b1af4e9 --- /dev/null +++ b/include/IRMosaicRefinementCommon.h @@ -0,0 +1,1719 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : mosaic_refinement_common.hxx +// Author : Pavel A. Koshevoy, Joel Spaltenstein +// Created : Mon Mar 26 10:34:39 MDT 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : Helper functions for automatic mosaic refinement. + +#ifndef MOSAIC_REFINEMENT_COMMON_HXX_ +#define MOSAIC_REFINEMENT_COMMON_HXX_ + +// the includes: +#include +#include +#include +#include +#include +#include + +#include "itkImageDuplicator.h" + +// system includes: +#include +#include +#include + +#ifndef WIN32 +# include +#endif + + +//---------------------------------------------------------------- +// regularize_displacements +// +extern void +regularize_displacements( // computed displacement vectors of the moving image + // grid transform control points, in mosaic space: + std::vector & xy_shift, + std::vector & mass, + + image_t::Pointer & dx, + image_t::Pointer & dy, + image_t::Pointer & db, + + // median filter radius: + const unsigned int & median_radius); + + +//---------------------------------------------------------------- +// refine_mosaic +// +template +void +refine_mosaic(the_log_t & log, + + // all the tiles in the mosaic, and their masks: + array2d(typename TImage::Pointer) & pyramid, + std::vector & mask, + std::vector & transform, + + unsigned int iterations_per_level) +{ + typedef itk::LinearInterpolateImageFunction interpolator_t; + + unsigned int pyramid_levels = pyramid.size(); + if (pyramid_levels == 0) + return; + + std::vector & image = pyramid[pyramid_levels - 1]; + unsigned int num_images = image.size(); + if (num_images == 0) + return; + + log << "iterations per level: " << iterations_per_level << endl + << "transform type: " << transform[0]->GetTransformTypeAsString() << endl; + + // try global refinement of the mosaic: + typedef itk::ImageMosaicVarianceMetric mosaic_metric_t; + + typename mosaic_metric_t::Pointer mosaic_metric = mosaic_metric_t::New(); + mosaic_metric->image_.resize(num_images); + mosaic_metric->mask_.resize(num_images); + mosaic_metric->transform_.resize(num_images); + for (unsigned int i = 0; i < num_images; i++) + { + mosaic_metric->image_[i] = image[i]; + mosaic_metric->mask_[i] = mask[i]; + mosaic_metric->transform_[i] = transform[i]; + } + + // FIXME: ITK doesn't have an API for this: + std::vector param_shared(transform[0]->GetNumberOfParameters(), false); + std::vector param_active(transform[0]->GetNumberOfParameters(), true); + + // setup the shared parameters mask: + mosaic_metric->setup_param_map(param_shared, param_active); + mosaic_metric->Initialize(); + + // setup the optimizer scales: + typename mosaic_metric_t::params_t parameter_scales = mosaic_metric->GetTransformParameters(); + parameter_scales.Fill(1.0); + + for (unsigned int level = 0; level < pyramid_levels && iterations_per_level > 0; level++) + { + for (unsigned int i = 0; i < num_images; i++) + { + mosaic_metric->image_[i] = pyramid[level][i]; + } + + typename mosaic_metric_t::measure_t metric_before = + mosaic_metric->GetValue(mosaic_metric->GetTransformParameters()); + + // run several iterations of the optimizer: + for (unsigned int k = 0; k < 3; k++) + { + typename mosaic_metric_t::params_t params_before = mosaic_metric->GetTransformParameters(); + + typename mosaic_metric_t::measure_t metric_after = + std::numeric_limits::max(); + + // use global refinement: + optimizer_t::Pointer optimizer = optimizer_t::New(); + typename optimizer_observer_t::Pointer observer = optimizer_observer_t::New(); + observer->log_ = &log; + optimizer->AddObserver(itk::IterationEvent(), observer); + optimizer->SetLog(&log); + optimizer->SetMinimize(true); + optimizer->SetNumberOfIterations(iterations_per_level); + optimizer->SetMinimumStepLength(1e-12); // min_step + optimizer->SetMaximumStepLength(1e-5); // max_step + optimizer->SetGradientMagnitudeTolerance(1e-6); + optimizer->SetRelaxationFactor(5e-1); + optimizer->SetCostFunction(mosaic_metric); + optimizer->SetInitialPosition(params_before); + optimizer->SetScales(parameter_scales); + optimizer->SetPickUpPaceSteps(5); + optimizer->SetBackTracking(true); + + // refine the mosaic: + try + { + log << "\n" << level << '.' << k << ": refining distortion transforms" << endl; + + optimizer->StartOptimization(); + } + catch (itk::ExceptionObject & exception) + { + // oops: + log << "optimizer threw an exception:" << endl << exception.what() << endl; + } + + mosaic_metric->SetTransformParameters(optimizer->GetBestParams()); + metric_after = optimizer->GetBestValue(); + + typename mosaic_metric_t::params_t params_after = mosaic_metric->GetTransformParameters(); + + log << "before: METRIC = " << metric_before << ", PARAMS = " << params_before << endl + << "after: METRIC = " << metric_after << ", PARAMS = " << params_after << endl; + + // quantify the improvement: + double improvement = 1.0 - metric_after / metric_before; + bool failed_to_improve = (metric_after - metric_before) >= 0.0; + bool negligible_improvement = !failed_to_improve && (improvement < 1e-3); + + if (!failed_to_improve) + { + log << "IMPROVEMENT: " << setw(3) << int(100.0 * improvement) << "%" << endl; + } + + if (failed_to_improve) + { + log << "NOTE: minimization failed, ignoring registration results..." << endl; + + // previous transform was better: + mosaic_metric->SetTransformParameters(params_before); + break; + } + else if (negligible_improvement) + { + log << "NOTE: improvement is negligible..." << endl; + break; + } + + // avoid recalculating the same metric: + metric_before = metric_after; + } + } +} + + +//---------------------------------------------------------------- +// calc_displacements +// +// This is the original version used for single-threaded execution +// +template +void +calc_displacements(the_log_t & log, + + // computed displacement vectors of the moving image + // grid transform control points, in mosaic space: + std::vector & xy_shift, + std::vector & mass, + + // a flag indicating whether the tiles + // have already been transformed into mosaic space: + bool tiles_already_warped, + + // fixed: + const TImage * tile_0, + const TMask * mask_0, + const base_transform_t * forward_0, + + // moving: + const TImage * tile_1, + const TMask * mask_1, + const itk::GridTransform * forward_1, + + // neighborhood size: + const unsigned int & neighborhood, + + // minimum acceptable neighborhood overlap ratio: + const double & min_overlap, + + // median filter radius: + const unsigned int & median_radius) +{ + // shortcuts: + const the_grid_transform_t & gt = forward_1->transform_; + unsigned int mesh_cols = gt.cols_ + 1; + unsigned int mesh_rows = gt.rows_ + 1; + unsigned mesh_size = gt.grid_.mesh_.size(); + xy_shift.assign(mesh_size, vec2d(0, 0)); + + // make sure both tiles have the same pixel spacing: + typename TImage::SpacingType sp = tile_1->GetSpacing(); + if (sp != tile_0->GetSpacing()) + return; + + // setup the local neighborhood: + typename TImage::SizeType sz; + sz[0] = neighborhood; + sz[1] = neighborhood; + + typename TImage::Pointer img[] = { make_image(sp, sz), make_image(sp, sz) }; + + typename TMask::Pointer msk[] = { make_image(sp, sz), make_image(sp, sz) }; + + // for each interpolation point, do a local neighborhood fft matching, + // and use the resulting displacement vector to adjust the mesh: + image_t::Pointer dx = make_image(mesh_cols, mesh_rows, 1.0, 0.0); + + image_t::Pointer dy = make_image(mesh_cols, mesh_rows, 1.0, 0.0); + + image_t::Pointer db = make_image(mesh_cols, mesh_rows, 1.0, 0.0); + + typename TImage::Pointer img_large; + typename TMask::Pointer msk_large; + + if (!tiles_already_warped) + { + typename TImage::SizeType sz_large(sz); + sz_large[0] *= 2; + sz_large[1] *= 2; + img_large = make_image(sp, sz_large); + msk_large = make_image(sp, sz_large); + } + + log << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl; +#pragma omp parallel for + for (int i = 0; i < (int)mesh_size; i++) + { + // shortcut: + const vertex_t & vertex = gt.grid_.mesh_[i]; + + // find the mosaic space coordinates of this vertex: + pnt2d_t center; + gt.transform_inv(vertex.uv_, center); + + // extract a neighborhood of the vertex from both tiles: + image_t::IndexType index; + index[0] = i % mesh_cols; + index[1] = i / mesh_cols; + dx->SetPixel(index, 0); + dy->SetPixel(index, 0); + db->SetPixel(index, 0); + + // feed the two neighborhoods into the FFT translation estimator: + vec2d_t shift(vec2d(0, 0)); + bool ok = tiles_already_warped ? + + refine_one_point_fft(log, + shift, + + tile_0, + mask_0, + + tile_1, + mask_1, + + center, + min_overlap, + + img[0].GetPointer(), + msk[0].GetPointer(), + + img[1].GetPointer(), + msk[1].GetPointer()) + : + + refine_one_point_fft(log, + shift, + + tile_0, + mask_0, + + tile_1, + mask_1, + + forward_0, + forward_1, + + center, + min_overlap, + + sz, + sp, + + img_large.GetPointer(), + msk_large.GetPointer(), + + img[0].GetPointer(), + msk[0].GetPointer(), + + img[1].GetPointer(), + msk[1].GetPointer()); + if (!ok) + { + continue; + } + + log << i << ". shift: " << shift << endl; + dx->SetPixel(index, shift[0]); + dy->SetPixel(index, shift[1]); + db->SetPixel(index, 1); + } + + // regularize the displacement vectors here: + regularize_displacements(xy_shift, mass, dx, dy, db, median_radius); +} + + +//---------------------------------------------------------------- +// calc_displacements_t +// +// +template +class calc_displacements_t : public the_transaction_t +{ +public: + calc_displacements_t( // a flag indicating whether the tiles + // have already been transformed into mosaic space: + bool tiles_already_warped, + + // fixed: + const TImage * tile_0, + const TMask * mask_0, + const base_transform_t * forward_0, + + // moving: + const TImage * tile_1, + const TMask * mask_1, + const itk::GridTransform * forward_1, + + // neighborhood size: + const unsigned int & neighborhood_size, + + // minimum acceptable neighborhood overlap ratio: + const double & min_overlap, + + // mesh node displacements: + image_t::Pointer dx, + image_t::Pointer dy, + + // mesh node displacement weights: + image_t::Pointer db, + + // mesh node index: + const std::list & index, + + // mesh node mosaic space coordinates: + const std::list & center) + : tiles_already_warped_(tiles_already_warped) + , + + tile_0_(tile_0) + , mask_0_(mask_0) + , forward_0_(forward_0) + , + + tile_1_(tile_1) + , mask_1_(mask_1) + , forward_1_(forward_1) + , + + min_overlap_(min_overlap) + , + + dx_(dx) + , dy_(dy) + , db_(db) + , + + index_(index.begin(), index.end()) + , center_(center.begin(), center.end()) + { + // make sure both tiles have the same pixel spacing: + sp_ = tile_1->GetSpacing(); + bool ok = (sp_ == tile_0->GetSpacing()); + if (!ok) + { + assert(false); + return; + } + + // setup the local neighborhood: + sz_[0] = neighborhood_size; + sz_[1] = neighborhood_size; + } + + // virtual: + void + execute(the_thread_interface_t * thread) + { + WRAP(the_terminator_t terminator("calc_displacements_t")); + + typename TImage::Pointer img[] = { make_image(sp_, sz_), make_image(sp_, sz_) }; + + typename TMask::Pointer msk[] = { make_image(sp_, sz_), make_image(sp_, sz_) }; + + typename TImage::Pointer img_large; + typename TMask::Pointer msk_large; + + if (!tiles_already_warped_) + { + typename TImage::SizeType sz_large(sz_); + sz_large[0] *= 2; + sz_large[1] *= 2; + img_large = make_image(sp_, sz_large); + msk_large = make_image(sp_, sz_large); + } + + std::size_t num_nodes = center_.size(); + for (std::size_t i = 0; i < num_nodes; i++) + { + // shortcuts: + const image_t::IndexType & index = index_[i]; + const pnt2d_t & center = center_[i]; + + // feed the two neighborhoods into the FFT translation estimator: + vec2d_t shift(vec2d(0, 0)); + + bool ok = tiles_already_warped_ ? + // if the tiles are already warped we don't + // need the transforms + refine_one_point_fft(*null_log(), + shift, + + tile_0_, + mask_0_, + + tile_1_, + mask_1_, + + center, + min_overlap_, + + img[0].GetPointer(), + msk[0].GetPointer(), + + img[1].GetPointer(), + msk[1].GetPointer()) + : + + refine_one_point_fft(*null_log(), + shift, + + tile_0_, + mask_0_, + + tile_1_, + mask_1_, + + forward_0_, + forward_1_, + + center, + min_overlap_, + + sz_, + sp_, + + img_large.GetPointer(), + msk_large.GetPointer(), + + img[0].GetPointer(), + msk[0].GetPointer(), + + img[1].GetPointer(), + msk[1].GetPointer()); + if (ok) + { + dx_->SetPixel(index, shift[0]); + dy_->SetPixel(index, shift[1]); + db_->SetPixel(index, 1); + } + } + } + + // a flag indicating whether the tiles + // have already been transformed into mosaic space: + bool tiles_already_warped_; + + // fixed: + const TImage * tile_0_; + const TMask * mask_0_; + const base_transform_t * forward_0_; + + // moving: + const TImage * tile_1_; + const TMask * mask_1_; + const itk::GridTransform * forward_1_; + + // moving tile spacing (must be the same as fixed tile spacing): + typename TImage::SpacingType sp_; + typename TImage::SizeType sz_; + + // minimum acceptable neighborhood overlap ratio: + const double min_overlap_; + + // mesh node displacements: + image_t::Pointer dx_; + image_t::Pointer dy_; + + // mesh node displacement weights: + image_t::Pointer db_; + + // mesh node index: + std::vector index_; + + // mesh node mosaic space coordinates: + std::vector center_; +}; + + +//---------------------------------------------------------------- +// calc_displacements_mt +// +// calculate transform mesh displacement vectors multi-threaded +// +template +void +calc_displacements_mt(unsigned int num_threads, + + the_log_t & log, + + // computed displacement vectors of the moving image + // grid transform control points, in mosaic space: + std::vector & xy_shift, + std::vector & mass, + + // a flag indicating whether the tiles + // have already been transformed into mosaic space: + bool tiles_already_warped, + + // fixed: + typename TImage::ConstPointer tile_0, + typename TMask::ConstPointer mask_0, + typename base_transform_t::ConstPointer forward_0, + + // moving: + typename TImage::ConstPointer tile_1, + typename TMask::ConstPointer mask_1, + typename itk::GridTransform::ConstPointer & forward_1, + + // neighborhood size: + const unsigned int & neighborhood_size, + + // minimum acceptable neighborhood overlap ratio: + const double & min_overlap, + + // median filter radius: + const unsigned int & median_radius) +{ + // make sure both tiles have the same pixel spacing: + if (tile_1->GetSpacing() != tile_0->GetSpacing()) + return; + + // shortcuts: + const the_grid_transform_t & gt = forward_1->transform_; + unsigned int mesh_cols = gt.cols_ + 1; + unsigned int mesh_rows = gt.rows_ + 1; + unsigned mesh_size = gt.grid_.mesh_.size(); + xy_shift.assign(mesh_size, vec2d(0, 0)); + + // for each interpolation point, do a local neighborhood fft matching, + // and use the resulting displacement vector to adjust the mesh: + image_t::Pointer dx = make_image(mesh_cols, mesh_rows, 1.0, 0.0); + + image_t::Pointer dy = make_image(mesh_cols, mesh_rows, 1.0, 0.0); + + image_t::Pointer db = make_image(mesh_cols, mesh_rows, 1.0, 0.0); + + the_thread_pool_t thread_pool(num_threads); + thread_pool.set_idle_sleep_duration(50); // 50 usec + + // split nodes between threads: + std::vector> node_index_list(num_threads); + std::vector> node_center_list(num_threads); + +#pragma omp parallel for + for (int i = 0; i < mesh_size; i++) + { + // shortcuts: + const vertex_t & vertex = gt.grid_.mesh_[i]; + const unsigned int which_thread = i % num_threads; + + // find the mosaic space coordinates of this vertex: + pnt2d_t center; + gt.transform_inv(vertex.uv_, center); + node_center_list[which_thread].push_back(center); + + // extract a neighborhood of the vertex from both tiles: + image_t::IndexType index; + index[0] = i % mesh_cols; + index[1] = i / mesh_cols; + dx->SetPixel(index, 0); + dy->SetPixel(index, 0); + db->SetPixel(index, 0); + + node_index_list[which_thread].push_back(index); + } + + // setup a transaction for each thread: + for (unsigned int i = 0; i < num_threads; i++) + { + calc_displacements_t * t = new calc_displacements_t(tiles_already_warped, + + tile_0, + mask_0, + forward_0, + + tile_1, + mask_1, + forward_1, + + neighborhood_size, + min_overlap, + + dx, + dy, + db, + + node_index_list[i], + node_center_list[i]); + + thread_pool.push_back(t); + } + + // run the transactions: + thread_pool.pre_distribute_work(); + suspend_itk_multithreading_t suspend_itk_mt; + thread_pool.start(); + thread_pool.wait(); + + // regularize the displacement vectors here: + regularize_displacements(xy_shift, mass, dx, dy, db, median_radius); +} + + +//---------------------------------------------------------------- +// intermediate_result_t +// +class intermediate_result_t +{ +public: + intermediate_result_t() {} + + intermediate_result_t(const intermediate_result_t & r) { *this = r; } + + intermediate_result_t(unsigned int num_neighbors, unsigned int mesh_rows, unsigned int mesh_cols) + : dx_(num_neighbors) + , dy_(num_neighbors) + , db_(num_neighbors) + { + for (unsigned int i = 0; i < num_neighbors; i++) + { + dx_[i] = make_image(mesh_cols, mesh_rows, 1.0, 0.0); + + dy_[i] = make_image(mesh_cols, mesh_rows, 1.0, 0.0); + + db_[i] = make_image(mesh_cols, mesh_rows, 1.0, 0.0); + } + } + + intermediate_result_t & + operator=(const intermediate_result_t & r) + { + if (this != &r) + { + dx_.resize(r.dx_.size()); + for (unsigned int i = 0; i < dx_.size(); i++) + { + dx_[i] = cast(r.dx_[i]); + } + + dy_.resize(r.dy_.size()); + for (unsigned int i = 0; i < dy_.size(); i++) + { + dy_[i] = cast(r.dy_[i]); + } + + db_.resize(r.db_.size()); + for (unsigned int i = 0; i < db_.size(); i++) + { + db_[i] = cast(r.db_[i]); + } + } + + return *this; + } + + std::vector dx_; + std::vector dy_; + std::vector db_; +}; + +//---------------------------------------------------------------- +// calc_intermediate_results_t +// +template +class calc_intermediate_results_t : public the_transaction_t +{ +public: + typedef typename img_ptr_t::ObjectType::Pointer::ObjectType TImage; + typedef typename msk_ptr_t::ObjectType::Pointer::ObjectType TMask; + + calc_intermediate_results_t(the_log_t & log, + unsigned int thread_offset, + unsigned int thread_stride, + std::vector & transform, + const std::vector & warped_tile, + const std::vector & warped_mask, + const std::vector> & neighbors, + const bool & tiles_already_warped, + const unsigned int & neighborhood_size, + const double & minimum_overlap, + const bool & keep_first_tile_fixed, + std::vector & results) + : + + log_(log) + , thread_offset_(thread_offset) + , thread_stride_(thread_stride) + , transform_(transform) + , warped_tile_(warped_tile) + , warped_mask_(warped_mask) + , neighbors_(neighbors) + , tiles_already_warped_(tiles_already_warped) + , neighborhood_(neighborhood_size) + , minimum_overlap_(minimum_overlap) + , keep_first_tile_fixed_(keep_first_tile_fixed) + , results_(results) + {} + + // virtual: + void + execute(the_thread_interface_t * thread) + { + WRAP(the_terminator_t terminator("calc_intermediate_results_t::execute")); + + unsigned int num_tiles = warped_tile_.size(); + unsigned int start = keep_first_tile_fixed_ ? 1 : 0; + + // setup the local neighborhood image buffers: + typename TImage::SpacingType sp = warped_tile_[start]->GetSpacing(); + typename TImage::SizeType sz; + sz[0] = neighborhood_; + sz[1] = neighborhood_; + + typename TImage::Pointer img[] = { make_image(sp, sz), make_image(sp, sz) }; + + typename TMask::Pointer msk[] = { make_image(sp, sz), make_image(sp, sz) }; + + typename TImage::Pointer img_large; + typename TMask::Pointer msk_large; + + if (!tiles_already_warped_) + { + typename TImage::SizeType sz_large(sz); + sz_large[0] *= 2; + sz_large[1] *= 2; + img_large = make_image(sp, sz_large); + msk_large = make_image(sp, sz_large); + } + + for (unsigned int tile_index = start; tile_index < num_tiles; tile_index++) + { + // shortcuts: + const unsigned int num_neighbors = neighbors_[tile_index].size(); + const TImage * tile = warped_tile_[tile_index]; + const TMask * mask = warped_mask_[tile_index]; + const itk::GridTransform * transform = transform_[tile_index]; + + const the_grid_transform_t & gt = transform->transform_; + const unsigned int mesh_cols = gt.cols_ + 1; + const unsigned int mesh_size = gt.grid_.mesh_.size(); + + // shortcuts to intermediate mesh refinement results for this tile: + intermediate_result_t & results = results_[tile_index]; + + for (unsigned int neighbor = 0; neighbor < num_neighbors; neighbor++) + { + WRAP(terminator.terminate_on_request()); + + // shortcuts: + image_t::Pointer & dx = results.dx_[neighbor]; + image_t::Pointer & dy = results.dy_[neighbor]; + image_t::Pointer & db = results.db_[neighbor]; + + const unsigned int neighbor_index = neighbors_[tile_index][neighbor]; + log_ << thread_offset_ << " thread, matching " << tile_index << ":" << neighbor_index << endl; + + // shortcuts: + const TImage * neighbor_tile = warped_tile_[neighbor_index]; + const TMask * neighbor_mask = warped_mask_[neighbor_index]; + const base_transform_t * neighbor_xform = transform_[neighbor_index]; + + for (unsigned int mesh_index = thread_offset_; mesh_index < mesh_size; mesh_index += thread_stride_) + { + // shortcut: + const vertex_t & vertex = gt.grid_.mesh_[mesh_index]; + + // find the mosaic space coordinates of this vertex: + pnt2d_t center; + gt.transform_inv(vertex.uv_, center); + + // extract a neighborhood of the vertex from both tiles: + image_t::IndexType index; + index[0] = mesh_index % mesh_cols; + index[1] = mesh_index / mesh_cols; + dx->SetPixel(index, 0); + dy->SetPixel(index, 0); + db->SetPixel(index, 0); + + // feed the two neighborhoods into the FFT translation estimator: + vec2d_t shift(vec2d(0, 0)); + bool ok = tiles_already_warped_ ? + + refine_one_point_fft(*null_log(), + shift, + + // fixed: + neighbor_tile, + neighbor_mask, + + // moving: + tile, + mask, + + center, + minimum_overlap_, + + img[0].GetPointer(), + msk[0].GetPointer(), + + img[1].GetPointer(), + msk[1].GetPointer()) + : + + refine_one_point_fft(*null_log(), + shift, + + // fixed: + neighbor_tile, + neighbor_mask, + + // moving: + tile, + mask, + + neighbor_xform, + transform, + + center, + minimum_overlap_, + + sz, + sp, + + img_large.GetPointer(), + msk_large.GetPointer(), + + img[0].GetPointer(), + msk[0].GetPointer(), + + img[1].GetPointer(), + msk[1].GetPointer()); + if (!ok) + { + continue; + } + + dx->SetPixel(index, shift[0]); + dy->SetPixel(index, shift[1]); + db->SetPixel(index, 1); + } + } + } + } + + the_log_t & log_; + unsigned int thread_offset_; + unsigned int thread_stride_; + std::vector & transform_; + const std::vector & warped_tile_; + const std::vector & warped_mask_; + const std::vector> & neighbors_; + const bool tiles_already_warped_; + const unsigned int neighborhood_; + const double minimum_overlap_; // neighborhood overlap + const bool keep_first_tile_fixed_; + std::vector & results_; +}; + + +//---------------------------------------------------------------- +// update_tile_mesh_t +// +class update_tile_mesh_t : public the_transaction_t +{ +public: + update_tile_mesh_t(the_log_t & log, + unsigned int tile_index, + bool keep_first_tile_fixed, + unsigned int median_filter_radius, + std::vector & transform, + std::vector & results) + : log_(log) + , tile_index_(tile_index) + , keep_first_tile_fixed_(keep_first_tile_fixed) + , median_filter_radius_(median_filter_radius) + , transform_(transform) + , results_(results) + {} + + // virtual: + void + execute(the_thread_interface_t * thread) + { + WRAP(the_terminator_t terminator("update_tile_mesh_t::execute")); + + log_ << tile_index_ << " mesh update" << endl; + + // shortcuts: + itk::GridTransform * transform = transform_[tile_index_]; + the_grid_transform_t & gt = transform->transform_; + const unsigned int mesh_size = gt.grid_.mesh_.size(); + + // shortcuts to intermediate mesh refinement results for this tile: + intermediate_result_t & results = results_[tile_index_]; + const unsigned int num_neighbors = results.dx_.size(); + + std::vector shift(mesh_size, vec2d(0, 0)); + std::vector mass(mesh_size, 0); + + for (unsigned int neighbor = 0; neighbor < num_neighbors; neighbor++) + { + WRAP(terminator.terminate_on_request()); + + // shortcuts: + image_t::Pointer & dx = results.dx_[neighbor]; + image_t::Pointer & dy = results.dy_[neighbor]; + image_t::Pointer & db = results.db_[neighbor]; + + std::vector neighbor_pull(mesh_size, vec2d(0, 0)); + regularize_displacements(neighbor_pull, mass, dx, dy, db, median_filter_radius_); + + // blend the displacement vectors: + for (unsigned int i = 0; i < mesh_size; i++) + { + shift[i] += neighbor_pull[i]; + } + } + + // FIXME: if (num_neighbors > 1) + if (!keep_first_tile_fixed_) + { + // normalize: + for (unsigned int i = 0; i < mesh_size; i++) + { + double scale = 1 / (1 + mass[i]); + shift[i] *= scale; + } + } + + gt.grid_.update(&(shift[0])); + transform->setup(gt); + } + + the_log_t & log_; + const unsigned int tile_index_; + const bool keep_first_tile_fixed_; + const unsigned int median_filter_radius_; + std::vector & transform_; + std::vector & results_; +}; + + +//---------------------------------------------------------------- +// refine_mosaic +// +// use FFT to refine the grid transforms directly: +// +template +void +refine_mosaic(the_log_t & log, + std::vector & transform, + const std::vector & tile, + const std::vector & mask, + const unsigned int & neighborhood, + const bool & prewarp_tiles = true, + const double & minimum_overlap = 0.25, + const unsigned int & median_radius = 1, + const unsigned int & num_passes = 1, + const bool & keep_first_tile_fixed = false) +{ + typedef itk::GridTransform itk_grid_transform_t; + + // shortcut: + unsigned int num_tiles = tile.size(); + if (num_tiles < 2) + return; + + std::vector ConstTile(tile.size()); + ConstTile.assign(tile.begin(), tile.end()); + + // image space bounding boxes: + std::vector image_min; + std::vector image_max; + calc_image_bboxes(ConstTile, image_min, image_max); + + std::vector ConstTransform(transform.size()); + ConstTransform.assign(transform.begin(), transform.end()); + + // mosaic space bounding boxes: + std::vector mosaic_min; + std::vector mosaic_max; + calc_mosaic_bboxes(ConstTransform, image_min, image_max, mosaic_min, mosaic_max, 16); + + unsigned int start = keep_first_tile_fixed ? 1 : 0; + + // find overlapping neighbors for each tile: + std::vector> neighbors(num_tiles); + for (unsigned int i = start; i < num_tiles; i++) + { + the_aa_bbox_t ibox; + ibox << p3x1_t(mosaic_min[i][0], mosaic_min[i][1], 0) << p3x1_t(mosaic_max[i][0], mosaic_max[i][1], 0); + + for (unsigned int j = 0; j < num_tiles; j++) + { + if (i == j) + continue; + + the_aa_bbox_t jbox; + jbox << p3x1_t(mosaic_min[j][0], mosaic_min[j][1], 0) << p3x1_t(mosaic_max[j][0], mosaic_max[j][1], 0); + + if (!ibox.intersects(jbox)) + continue; + + neighbors[i].push_back(j); + } + } + + typedef itk::ImageDuplicator image_duplicator_t; + typedef itk::ImageDuplicator mask_duplicator_t; + + std::vector warped_tile(num_tiles); + std::vector warped_mask(num_tiles); + + if (keep_first_tile_fixed) + { + typename image_duplicator_t::Pointer image_duplicator = image_duplicator_t::New(); + typename mask_duplicator_t::Pointer mask_duplicator = mask_duplicator_t::New(); + + image_duplicator->SetInputImage(tile[0]); + mask_duplicator->SetInputImage(mask[0]); + + image_duplicator->Update(); + mask_duplicator->Update(); + + warped_tile[0] = image_duplicator->GetOutput(); + warped_mask[0] = mask_duplicator->GetOutput(); + } + + for (unsigned int pass = 0; pass < num_passes; pass++) + { + log << "--------------------------- pass " << pass << " ---------------------------" << endl; + + if (prewarp_tiles) + { + // warp the tiles: + for (unsigned int i = start; i < num_tiles; i++) + { + log << setw(4) << i << ". warping image tile" << endl; + warped_tile[i] = warp((typename image_t::ConstPointer)tile[i], transform[i].GetPointer()); + + if (mask[i].GetPointer() != NULL) + { + log << " warping image tile mask" << endl; + warped_mask[i] = warp((typename mask_t::ConstPointer)mask[i], transform[i].GetPointer()); + } + } + } + + std::vector> shift(num_tiles); + for (unsigned int i = start; i < num_tiles; i++) + { + the_grid_transform_t & gt = transform[i]->transform_; + unsigned int mesh_size = gt.grid_.mesh_.size(); + + std::vector> shift_i(neighbors[i].size()); + std::vector mass(mesh_size, 0); + + for (unsigned int k = 0; k < neighbors[i].size(); k++) + { + unsigned int j = neighbors[i][k]; + log << "matching " << i << ":" << j << endl; + + calc_displacements(*null_log(), + shift_i[k], + mass, + + // tiles-already-warped flag: + prewarp_tiles, + + // fixed: + warped_tile[j], + warped_mask[j], + transform[j].GetPointer(), + + // moving: + warped_tile[i], + warped_mask[i], + transform[i].GetPointer(), + + neighborhood, + minimum_overlap, + median_radius); + } + + // blend the displacement vectors: + shift[i].assign(mesh_size, vec2d(0, 0)); + + for (unsigned int j = 0; j < shift_i.size(); j++) + { + for (unsigned int k = 0; k < mesh_size; k++) + { + shift[i][k] += shift_i[j][k]; + } + } + + if (!keep_first_tile_fixed) + { + // normalize: + for (unsigned int k = 0; k < mesh_size; k++) + { + double scale = 1 / (1 + mass[k]); + shift[i][k] *= scale; + } + } + } + + for (unsigned int i = start; i < num_tiles; i++) + { + the_grid_transform_t & gt = transform[i]->transform_; + gt.grid_.update(&(shift[i][0])); + transform[i]->setup(gt); + } + } +} + + +//---------------------------------------------------------------- +// warp_tile_transaction_t +// +template +class warp_tile_transaction_t : public the_transaction_t +{ +public: + warp_tile_transaction_t(the_log_t & log, + unsigned int tile_index, + itk::GridTransform::Pointer & transform, + const typename image_t::ConstPointer & tile, + const typename mask_t::ConstPointer & mask, + std::vector & warped_tile, + std::vector & warped_mask) + : log_(log) + , tile_index_(tile_index) + , transform_(transform) + , tile_(tile) + , mask_(mask) + , warped_tile_(warped_tile) + , warped_mask_(warped_mask) + {} + + // virtual: + void + execute(the_thread_interface_t * thread) + { + WRAP(the_terminator_t terminator("warp_tile_transaction_t")); + + log_ << setw(4) << tile_index_ << ". warping image tile" << endl; + warped_tile_[tile_index_] = warp(tile_, transform_.GetPointer()); + + if (mask_.GetPointer() != NULL) + { + log_ << setw(4) << tile_index_ << ". warping image tile mask" << endl; + warped_mask_[tile_index_] = warp(mask_, transform_.GetPointer()); + } + } + + the_log_t & log_; + unsigned int tile_index_; + itk::GridTransform::Pointer & transform_; + const typename image_t::ConstPointer & tile_; + const typename mask_t::ConstPointer & mask_; + std::vector & warped_tile_; + std::vector & warped_mask_; +}; + + +//---------------------------------------------------------------- +// refine_one_tile_t +// +template +class refine_one_tile_t : public the_transaction_t +{ +public: + refine_one_tile_t(the_log_t & log, + unsigned int tile_index, + std::vector & transform, + const std::vector & warped_tile, + const std::vector & warped_mask, + const std::vector> & neighbors, + const bool & tiles_already_warped, + const unsigned int & neighborhood_size, + const double & minimum_overlap, // neighbrhood overlap + const unsigned int & median_filter_radius, // for outliers + const bool & keep_first_tile_fixed, + std::vector> & shift) + : + + log_(log) + , tile_index_(tile_index) + , transform_(transform) + , warped_tile_(warped_tile) + , warped_mask_(warped_mask) + , neighbors_(neighbors) + , tiles_already_warped_(tiles_already_warped) + , neighborhood_(neighborhood_size) + , minimum_overlap_(minimum_overlap) + , median_radius_(median_filter_radius) + , keep_first_tile_fixed_(keep_first_tile_fixed) + , shift_(shift) + {} + + // virtual: + void + execute(the_thread_interface_t * thread) + { + WRAP(the_terminator_t terminator("refine_one_tile_t")); + + the_grid_transform_t & gt = transform_[tile_index_]->transform_; + unsigned int mesh_size = gt.grid_.mesh_.size(); + std::vector mass(mesh_size, 0); + + unsigned int num_neighbors = neighbors_[tile_index_].size(); + std::vector> shift_i(num_neighbors); + + for (unsigned int k = 0; k < num_neighbors; k++) + { + unsigned int j = neighbors_[tile_index_][k]; + log_ << "matching " << tile_index_ << ":" << j << endl; + + calc_displacements(*null_log(), + shift_i[k], + mass, + + // tiles-already-warped flag: + tiles_already_warped_, + + // fixed: + warped_tile_[j], + warped_mask_[j], + transform_[j].GetPointer(), + + // moving: + warped_tile_[tile_index_], + warped_mask_[tile_index_], + transform_[tile_index_].GetPointer(), + + neighborhood_, + minimum_overlap_, + median_radius_); + } + + // blend the displacement vectors: + shift_[tile_index_].assign(mesh_size, vec2d(0, 0)); + + for (unsigned int j = 0; j < num_neighbors; j++) + { + for (unsigned int k = 0; k < mesh_size; k++) + { + shift_[tile_index_][k] += shift_i[j][k]; + } + } + + if (!keep_first_tile_fixed_) + { + // normalize: + for (unsigned int k = 0; k < mesh_size; k++) + { + double scale = 1 / (1 + mass[k]); + shift_[tile_index_][k] *= scale; + } + } + } + + the_log_t & log_; + unsigned int tile_index_; + std::vector & transform_; + const std::vector & warped_tile_; + const std::vector & warped_mask_; + const std::vector> & neighbors_; + const bool tiles_already_warped_; + const unsigned int neighborhood_; + const double minimum_overlap_; // neighbrhood overlap + const unsigned int median_radius_; // for outliers + const bool keep_first_tile_fixed_; + std::vector> & shift_; +}; + + +//---------------------------------------------------------------- +// refine_mosaic_mt +// +template +void +refine_mosaic_mt(the_log_t & log, // text output stream + std::vector & transform, + const std::vector & tile, + const std::vector & mask, + const unsigned int & neighborhood_size, + const bool & prewarp_tiles, // FIXME: always true? + const double & minimum_overlap, // neighbrhood overlap + const unsigned int & median_radius, // for outliers + const unsigned int & num_passes, + const bool & keep_first_tile_fixed, // FIXME: stos only? + const double & displacement_threshold, + unsigned int num_threads) // max concurrent threads +{ + if (num_threads == 1) + { + refine_mosaic(log, + transform, + tile, + mask, + neighborhood_size, + prewarp_tiles, + minimum_overlap, + median_radius, + num_passes, + keep_first_tile_fixed); + return; + } + + typedef itk::GridTransform itk_grid_transform_t; + + // shortcut: + unsigned int num_tiles = tile.size(); + if (num_tiles < 2) + return; + + log << "num tiles: " << num_tiles; + // image space bounding boxes: + std::vector ConstTile; + std::vector image_min; + std::vector image_max; + + ConstTile.reserve(tile.size()); + ConstTile.assign(tile.begin(), tile.end()); + + std::vector ConstTransforms; + ConstTransforms.reserve(transform.size()); + + ConstTransforms.assign(transform.begin(), transform.end()); + + calc_image_bboxes(ConstTile, image_min, image_max); + + // mosaic space bounding boxes: + std::vector mosaic_min; + std::vector mosaic_max; + calc_mosaic_bboxes(ConstTransforms, image_min, image_max, mosaic_min, mosaic_max, 16); + + // Relative to the image size. + // double threshold = + // displacement_threshold * std::abs(mosaic_min[0][0] - mosaic_max[0][0]); + + // Relative to a single pixel. + double threshold = displacement_threshold; + + unsigned int start = keep_first_tile_fixed ? 1 : 0; + + // find overlapping neighbors for each tile: + std::vector> neighbors(num_tiles); + for (unsigned int i = start; i < num_tiles; i++) + { + the_aa_bbox_t ibox; + ibox << p3x1_t(mosaic_min[i][0], mosaic_min[i][1], 0) << p3x1_t(mosaic_max[i][0], mosaic_max[i][1], 0); + + for (unsigned int j = 0; j < num_tiles; j++) + { + if (i == j) + continue; + + the_aa_bbox_t jbox; + jbox << p3x1_t(mosaic_min[j][0], mosaic_min[j][1], 0) << p3x1_t(mosaic_max[j][0], mosaic_max[j][1], 0); + + if (!ibox.intersects(jbox)) + continue; + + neighbors[i].push_back(j); + } + } + + std::vector warped_tile(num_tiles); + std::vector warped_mask(num_tiles); + + double last_average = std::numeric_limits::max(); + + typedef itk::ImageDuplicator image_duplicator_t; + typedef itk::ImageDuplicator mask_duplicator_t; + + // Initialize "warped" tiles. + // If warping is actually requested do it on each pass + for (int i = 0; i < num_tiles; i++) + { + typename image_duplicator_t::Pointer image_duplicator = image_duplicator_t::New(); + typename mask_duplicator_t::Pointer mask_duplicator = mask_duplicator_t::New(); + + image_duplicator->SetInputImage(tile[i]); + image_duplicator->Update(); + warped_tile[i] = image_duplicator->GetOutput(); + + if (mask[i].GetPointer() != NULL) + { + mask_duplicator->SetInputImage(mask[i]); + mask_duplicator->Update(); + warped_mask[i] = mask_duplicator->GetOutput(); + } + } + + the_thread_pool_t thread_pool(num_threads); + thread_pool.set_idle_sleep_duration(50); // 50 usec + + for (unsigned int pass = 0; pass < num_passes; pass++) + { + double major_percent = 0.15 + 0.8 * ((double)pass / (double)num_passes); + double next_major = 0.15 + 0.8 * ((double)(pass + 1) / (double)num_passes); + set_major_progress(major_percent); + + log << "--------------------------- pass " << pass << " ---------------------------" << endl; + + if (prewarp_tiles) + { + // warp the tiles -- split into a set of transactions and executed: + + std::list schedule; + for (unsigned int i = start; i < num_tiles; i++) + { + warp_tile_transaction_t * t = new warp_tile_transaction_t( + log, i, transform[i], tile[i], mask[i], warped_tile, warped_mask); + schedule.push_back(t); + } + + thread_pool.push_back(schedule); + thread_pool.pre_distribute_work(); + suspend_itk_multithreading_t suspend_itk_mt; + thread_pool.start(); + thread_pool.wait(); + } + + set_minor_progress(0.2, next_major); + +#if 1 + // calculating displacements: + std::vector> shift(num_tiles); + + // this is the original coarse-scale parallelization: + unsigned int num_tiles_distributed = (num_tiles - start) - (num_tiles - start) % num_threads; + + unsigned int num_tiles_remaining = (num_tiles - start) - num_tiles_distributed; + + std::list schedule; + for (unsigned int i = 0; i < num_tiles_distributed; i++) + { + unsigned int index = start + i; + + typedef refine_one_tile_t tile_transaction_t; + + tile_transaction_t * t = new tile_transaction_t(log, + index, + transform, + warped_tile, + warped_mask, + neighbors, + prewarp_tiles, + neighborhood_size, + minimum_overlap, + median_radius, + keep_first_tile_fixed, + shift); + schedule.push_back(t); + } + + thread_pool.push_back(schedule); + // thread_pool.pre_distribute_work(); + suspend_itk_multithreading_t suspend_itk_mt; + thread_pool.start(); + thread_pool.wait(); + + set_minor_progress(0.9, next_major); + + // this is the broken fine-scale parallelization: + for (unsigned int i = 0; i < num_tiles_remaining; i++) + { + unsigned int index = start + num_tiles_distributed + i; + + the_grid_transform_t & gt = transform[index]->transform_; + unsigned int mesh_size = gt.grid_.mesh_.size(); + + std::vector> shift_i(neighbors[index].size()); + std::vector mass(mesh_size, 0); + + for (unsigned int k = 0; k < neighbors[index].size(); k++) + { + unsigned int j = neighbors[index][k]; + log << "matching " << index << ":" << j << endl; + + calc_displacements(*null_log(), + shift_i[k], + mass, + + // tiles-already-warped flag: + prewarp_tiles, + + // fixed: + (typename image_t::ConstPointer)warped_tile[j], + (typename mask_t::ConstPointer)warped_mask[j], + (base_transform_t::ConstPointer)transform[j].GetPointer(), + + // moving: + (typename image_t::ConstPointer)warped_tile[index], + (typename mask_t::ConstPointer)warped_mask[index], + (itk::GridTransform::ConstPointer)transform[index].GetPointer(), + + neighborhood_size, + minimum_overlap, + median_radius); + } + + // blend the displacement vectors: + shift[index].assign(mesh_size, vec2d(0, 0)); + + for (unsigned int j = 0; j < shift_i.size(); j++) + { + for (unsigned int k = 0; k < mesh_size; k++) + { + shift[index][k] += shift_i[j][k]; + } + } + + if (!keep_first_tile_fixed) + { + // normalize: + for (unsigned int k = 0; k < mesh_size; k++) + { + double scale = 1 / (1 + mass[k]); + shift[index][k] *= scale; + } + } + } + + // once all displacement calculations are + // finished the transform grids can be updated and we + // can move on to the next pass: + + for (unsigned int i = start; i < num_tiles; i++) + { + the_grid_transform_t & gt = transform[i]->transform_; + gt.grid_.update(&(shift[i][0])); + transform[i]->setup(gt); + } + +#else + // this is the "improved" fine scale parallelization: + + // setup intermediate mesh refinement result structures: + std::vector results(num_tiles); + for (unsigned int i = start; i < num_tiles; i++) + { + const unsigned int num_neighbors = neighbors[i].size(); + const the_grid_transform_t & gt = transform[i]->transform_; + const unsigned int mesh_rows = gt.rows_ + 1; + const unsigned int mesh_cols = gt.cols_ + 1; + + intermediate_result_t & result = results[i]; + result = intermediate_result_t(num_neighbors, mesh_rows, mesh_cols); + } + + std::list schedule; + for (unsigned int i = 0; i < num_threads; i++) + { + typedef calc_intermediate_results_t tile_transaction_t; + + tile_transaction_t * t = new tile_transaction_t(log, + i, + num_threads, + transform, + warped_tile, + warped_mask, + neighbors, + prewarp_tiles, + neighborhood_size, + minimum_overlap, + keep_first_tile_fixed, + results); + schedule.push_back(t); + } + + thread_pool.push_back(schedule); + thread_pool.pre_distribute_work(); + suspend_itk_multithreading_t suspend_itk_mt; + thread_pool.start(); + thread_pool.wait(); + + for (unsigned int i = start; i < num_tiles; i++) + { + update_tile_mesh_t * t = new update_tile_mesh_t(log, i, keep_first_tile_fixed, median_radius, transform, results); + schedule.push_back(t); + } + + thread_pool.push_back(schedule); + thread_pool.pre_distribute_work(); + thread_pool.start(); + thread_pool.wait(); +#endif + double worst = 0.0, avg = 0.0, count = 0.0; + for (int i = 0; i < (int)shift.size(); i++) + { + for (int k = 0; k < (int)shift[i].size(); k++) + { + if (std::abs(shift[i][k][0]) > worst) + worst = std::abs(shift[i][k][0]); + if (std::abs(shift[i][k][1]) > worst) + worst = std::abs(shift[i][k][1]); + + avg += std::abs(shift[i][k][0]) + std::abs(shift[i][k][1]); + count += 2; + } + } + avg /= count; + cout << pass << " Average Displacement: " << avg << " Max Displacement: " << worst << endl; + + // If there's an exact cutoff... + if (count > 0) + { + if (avg <= threshold) + break; + else if (avg >= last_average) + break; + last_average = avg; + } + } +} + + +#endif // MOSAIC_REFINEMENT_COMMON_HXX_ diff --git a/include/IROptimizerObserver.h b/include/IROptimizerObserver.h new file mode 100644 index 0000000..d78bd90 --- /dev/null +++ b/include/IROptimizerObserver.h @@ -0,0 +1,74 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : optimizer_observer.hxx +// Author : Pavel A. Koshevoy +// Created : Fri Mar 23 12:47:03 MDT 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An optimizer observer class. + +#ifndef OPTIMIZER_OBSERVER_HXX_ +#define OPTIMIZER_OBSERVER_HXX_ + +// the includes: +#include "utils/the_log.hxx" + +// ITK includes: +#include + + +//---------------------------------------------------------------- +// optimizer_observer_t +// +template +class optimizer_observer_t : public itk::Command +{ +public: + typedef optimizer_observer_t Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer Pointer; + + itkNewMacro(Self); + + void Execute(itk::Object *caller, const itk::EventObject & event) + { Execute((const itk::Object *)(caller), event); } + + void Execute(const itk::Object * object, const itk::EventObject & event) + { + if (typeid(event) != typeid(itk::IterationEvent)) return; + + const TOptimizer * optimizer = dynamic_cast(object); + *log_ << (unsigned int)(optimizer->GetCurrentIteration()) << '\t' + << optimizer->GetValue() << '\t' + << optimizer->GetCurrentPosition() << endl; + } + + the_log_t * log_; + +protected: + optimizer_observer_t(): + log_(cerr_log()) + {} +}; + + +#endif // OPTIMIZER_OBSERVER_HXX_ diff --git a/include/itkGridTransform.h b/include/itkGridTransform.h new file mode 100644 index 0000000..d912f55 --- /dev/null +++ b/include/itkGridTransform.h @@ -0,0 +1,325 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkGridTransform.h +// Author : Pavel A. Koshevoy +// Created : 2006/11/20 22:48 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A discontinuous transform -- a uniform grid of vertices is +// mapped to an image. At each vertex, in addition to image +// space coordinates, a second set of coordinates is stored. +// This is similar to texture mapped OpenGL triangle meshes, +// where the texture coordinates correspond to the image space +// vertex coordinates. + +#ifndef __itkGridTransform_h +#define __itkGridTransform_h + +// system includes: +#include +#include +#include + +// ITK includes: +#include +#include +#include +#include +#include + +// local includes: +#include "itk/the_grid_transform.hxx" + + +//---------------------------------------------------------------- +// itk::GridTransform +// +namespace itk +{ + class GridTransform : public Transform + { + public: + // standard typedefs: + typedef GridTransform Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + typedef Transform Superclass; + + // Base inverse transform type: + typedef Superclass::InverseTransformType InverseTransformType; + typedef SmartPointer< InverseTransformType > InverseTransformPointer; + + // RTTI: + itkTypeMacro(GridTransform, Transform); + + // macro for instantiation through the object factory: + itkNewMacro(Self); + + /** Standard scalar type for this class. */ + typedef double ScalarType; + + /** Dimension of the domain space. */ + itkStaticConstMacro(InputSpaceDimension, unsigned int, 2); + itkStaticConstMacro(OutputSpaceDimension, unsigned int, 2); + + // shortcuts: + typedef Superclass::ParametersType ParametersType; + typedef Superclass::JacobianType JacobianType; + + typedef Superclass::InputPointType InputPointType; + typedef Superclass::OutputPointType OutputPointType; + + // virtual: + OutputPointType TransformPoint(const InputPointType & x) const + { + OutputPointType y; + if (is_inverse()) + { + pnt2d_t uv; + uv[0] = (x[0] - transform_.tile_min_[0]) / transform_.tile_ext_[0]; + uv[1] = (x[1] - transform_.tile_min_[1]) / transform_.tile_ext_[1]; + transform_.transform_inv(uv, y); + } + else + { + transform_.transform(x, y); + y[0] *= transform_.tile_ext_[0]; + y[1] *= transform_.tile_ext_[1]; + y[0] += transform_.tile_min_[0]; + y[1] += transform_.tile_min_[1]; + } + + // ITK does not handle NaN numbers: + if (y[0] != y[0]) + { + y[0] = std::numeric_limits::max(); + y[1] = y[0]; + } + + return y; + } + + // Inverse transformations: + // If y = Transform(x), then x = BackTransform(y); + // if no mapping from y to x exists, then an exception is thrown. + InputPointType BackTransformPoint(const OutputPointType & y) const + { + InputPointType x; + if (is_inverse()) + { + transform_.transform(y, x); + x[0] *= transform_.tile_ext_[0]; + x[1] *= transform_.tile_ext_[1]; + x[0] += transform_.tile_min_[0]; + x[1] += transform_.tile_min_[1]; + } + else + { + pnt2d_t uv; + uv[0] = (y[0] - transform_.tile_min_[0]) / transform_.tile_ext_[0]; + uv[1] = (y[1] - transform_.tile_min_[1]) / transform_.tile_ext_[1]; + transform_.transform_inv(uv, x); + } + + // ITK does not handle NaN numbers: + if (x[0] != x[0]) + { + x[0] = std::numeric_limits::max(); + x[1] = x[0]; + } + + return x; + } + + // virtual: + void SetFixedParameters(const ParametersType & params) + { this->m_FixedParameters = params; } + + // virtual: + const ParametersType & GetFixedParameters() const + { + ParametersType params = this->m_FixedParameters; + + params[1] = transform_.rows_; + params[2] = transform_.cols_; + params[3] = transform_.tile_min_[0]; + params[4] = transform_.tile_min_[1]; + params[5] = transform_.tile_ext_[0]; + params[6] = transform_.tile_ext_[1]; + + // update the parameters vector: + GridTransform * fake = const_cast(this); + fake->m_FixedParameters = params; + return this->m_FixedParameters; + } + + // virtual: + void SetParameters(const ParametersType & params) + { + this->m_Parameters = params; + + size_t rows = int(this->m_FixedParameters[1]); + size_t cols = int(this->m_FixedParameters[2]); + + pnt2d_t tile_min; + tile_min[0] = this->m_FixedParameters[3]; + tile_min[1] = this->m_FixedParameters[4]; + + pnt2d_t tile_max; + tile_max[0] = tile_min[0] + this->m_FixedParameters[5]; + tile_max[1] = tile_min[1] + this->m_FixedParameters[6]; + + std::vector xy((rows + 1) * (cols + 1)); + + unsigned int num_points = xy.size(); + assert(2 * num_points == params.size()); + for (unsigned int i = 0; i < num_points; i++) + { + xy[i][0] = params[i * 2]; + xy[i][1] = params[i * 2 + 1]; + } + + transform_.setup(rows, + cols, + tile_min, + tile_max, + xy); + } + + // virtual: + const ParametersType & GetParameters() const + { + ParametersType params(GetNumberOfParameters()); + unsigned int num_pts = params.size() / 2; + unsigned int num_cols = (transform_.cols_ + 1); + for (unsigned int i = 0; i < num_pts; i++) + { + unsigned int row = i / num_cols; + unsigned int col = i % num_cols; + + const pnt2d_t & xy = transform_.vertex(row, col).xy_; + unsigned int idx = 2 * i; + params[idx] = xy[0]; + params[idx + 1] = xy[1]; + } + + // update the parameters vector: + GridTransform * fake = const_cast(this); + fake->m_Parameters = params; + return this->m_Parameters; + } + + // virtual: + unsigned int GetNumberOfParameters() const + { return 2 * transform_.grid_.mesh_.size(); } + + // virtual: return an inverse of this transform. + InverseTransformPointer GetInverse() const + { + GridTransform::Pointer inv = GridTransform::New(); + inv->setup(transform_, !is_inverse()); + return inv.GetPointer(); + } + + // setup the transform: + void setup(const the_grid_transform_t & transform, + const bool & is_inverse = false) + { + transform_ = transform; + GetParameters(); + GetFixedParameters(); + this->m_Jacobian.SetSize(2, GetNumberOfParameters()); + this->m_FixedParameters[0] = is_inverse ? 1.0 : 0.0; + } + + // inverse transform flag check: + inline bool is_inverse() const + { return this->m_FixedParameters[0] != 0.0; } + + // virtual: + const JacobianType & + GetJacobian(const InputPointType & point) const + { + // FIXME: 20061227 -- this function was written and not tested: + + // these scales are necessary to account for the fact that + // the_grid_transform_t expects uv in the [0,1]x[0,1] range, + // where as we remap it into the image tile physical coordinates + // according to tile_min_ and tile_ext_: + double Su = transform_.tile_ext_[0]; + double Sv = transform_.tile_ext_[1]; + + unsigned int idx[3]; + double jac[12]; + this->m_Jacobian.SetSize(2, GetNumberOfParameters()); + this->m_Jacobian.Fill(0.0); + if (transform_.jacobian(point, idx, jac)) + { + for (unsigned int i = 0; i < 3; i++) + { + unsigned int addr = idx[i] * 2; + this->m_Jacobian(0, addr) = Su * jac[i * 2]; + this->m_Jacobian(0, addr + 1) = Su * jac[i * 2 + 1]; + this->m_Jacobian(1, addr) = Sv * jac[i * 2 + 6]; + this->m_Jacobian(1, addr + 1) = Sv * jac[i * 2 + 7]; + } + } + + return this->m_Jacobian; + } + + protected: + GridTransform(): + Superclass(2, 0) + { + this->m_FixedParameters.SetSize(7); + + // initialize the inverse flag: + this->m_FixedParameters[0] = 0.0; + + // grid size: + this->m_FixedParameters[1] = 0.0; + this->m_FixedParameters[2] = 0.0; + + // tile bbox: + this->m_FixedParameters[3] = std::numeric_limits::max(); + this->m_FixedParameters[4] = this->m_FixedParameters[3]; + this->m_FixedParameters[5] = -(this->m_FixedParameters[3]); + this->m_FixedParameters[6] = -(this->m_FixedParameters[3]); + } + + private: + // disable default copy constructor and assignment operator: + GridTransform(const Self & other); + const Self & operator = (const Self & t); + + public: + // the actual transform: + the_grid_transform_t transform_; + + }; // class GridTransform + +} // namespace itk + + +#endif // __itkGridTransform_h diff --git a/include/itkImageMosaicVarianceMetric.h b/include/itkImageMosaicVarianceMetric.h new file mode 100644 index 0000000..d5d9213 --- /dev/null +++ b/include/itkImageMosaicVarianceMetric.h @@ -0,0 +1,317 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkImageMosaicVarianceMetric.h +// Author : Pavel A. Koshevoy +// Created : 2005/06/22 17:19 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A metric class measuring mean pixel variance within +// the mosaic tile overlap regions. The metric derivative + +#ifndef __itkImageMosaicVarianceMetric_h +#define __itkImageMosaicVarianceMetric_h + +#include +#include +#include +#include +#include +#include + +/** . + This class computes the mean pixel variance across several images + within the overlapping regions of the mosaic. Each image is warped + by a corresponding transform. All of the transforms must be of the + same type (and therefore have the same number of parameters). Some + of the transforms parameters may be shared across transforms. The + shared/unique parameters are specified by a bit vector + (true - shared, false - unique). This is usefull for radial + distortion transforms where the translation parameters are unique + for each image but the distortion parameters are the same across + all images (all images are assumed to have been distorted by + the same lens). + */ + +namespace itk +{ + +/** \class ImageMosaicVarianceMetric + * \brief Computes mean pixel variance within the overlapping regions + * of a mosaic. + * + */ +template +class ITK_EXPORT ImageMosaicVarianceMetric : public SingleValuedCostFunction +{ +public: + /** Standard class typedefs. */ + typedef ImageMosaicVarianceMetric Self; + typedef SingleValuedCostFunction Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(ImageMosaicVarianceMetric, SingleValuedCostFunction); + + /** Type of the moving Image. */ + typedef TInterpolator interpolator_t; + typedef TImage image_t; + typedef typename TImage::PixelType pixel_t; + + /** Constant for the image dimension */ + itkStaticConstMacro(ImageDimension, unsigned int, TImage::ImageDimension); + + typedef itk::Image mask_t; + + /** Type of the transform base class */ + typedef Transform + transform_t; + + typedef typename transform_t::InputPointType point_t; + typedef typename transform_t::ParametersType params_t; + typedef typename transform_t::JacobianType jacobian_t; + + /** Gaussian filter to compute the gradient of the mosaic images */ + // typedef typename NumericTraits::RealType + typedef float scalar_t; + + typedef CovariantVector gradient_pixel_t; + + typedef Image gradient_image_t; + + typedef GradientRecursiveGaussianImageFilter gradient_filter_t; + + /** Type of the measure. */ + typedef Superclass::MeasureType measure_t; + + /** Type of the derivative. */ + typedef Superclass::DerivativeType derivative_t; + + /** Get the number of pixels considered in the computation. */ + itkGetConstReferenceMacro(NumberOfPixelsCounted, unsigned long); + + /** Set/Get the region over which the metric will be computed */ + itkSetMacro(MosaicRegion, typename image_t::RegionType); + itkGetConstReferenceMacro(MosaicRegion, typename image_t::RegionType); + + /** Initialize transform parameter offsets. + + Setup the mapping from individual transform parameter indices to + the corresponding indices in the concatenated parameter vector. + + The radial distortion transforms typically share the same parameters + across all images in the mosaic, therefore it is slow and unnecessary + to duplicate the same parameters parameters in the concatenated + parameters vector. Instead, the parameters are store only once. + This also allows a significant speedup of the optimization. On the + other hand, translation parameters are typically unique for each + image. + + params_shared is a vector of boolean flags indicating which of + the transform parameters are shared (true) or unique (false). + + NOTE: this function will fail if the transform_ has not been + initialized. + */ + void + setup_param_map(const std::vector & params_shared, const std::vector & params_active); + + /** Concatenate parameters of each transform into one big vector. + + NOTE: this function will fail if SetupTransformParameterOffsets + has not been called yet. + */ + params_t + GetTransformParameters() const; + + /** virtual: Set the parameters defining the transforms: + + NOTE: this function will fail if SetupTransformParameterOffsets + has not been called yet. + */ + void + SetTransformParameters(const params_t & parameters) const; + + /** virtual: Return the number of parameters required by the transforms. + + NOTE: this function will fail if SetupTransformParameterOffsets + has not been called yet. + */ + unsigned int + GetNumberOfParameters() const; + + /** Initialize the Metric by making sure that all the components + are present and plugged together correctly. + */ + virtual void + Initialize() throw(ExceptionObject); + + /** virtual: Get the value for single valued optimizers. */ + measure_t + GetValue(const params_t & parameters) const + { + measure_t measure; + derivative_t derivative; + GetValueAndDerivative(parameters, measure, derivative); + return measure; + } + + /** virtual: Get the derivatives of the match measure. */ + void + GetDerivative(const params_t & parameters, derivative_t & derivative) const + { + measure_t measure; + GetValueAndDerivative(parameters, measure, derivative); + } + + /** virtual: Get value and derivatives for multiple valued optimizers. */ + void + GetValueAndDerivative(const params_t & parameters, measure_t & value, derivative_t & derivative) const; + + //---------------------------------------------------------------- + // image_data_t + // + // This datastructure is used to speed up access to relevant + // information when evaluating GetValueAndDerivative: + // + class image_data_t + { + public: + image_data_t() + : id_(~0) + , P_(measure_t(0)) + , dPdx_(measure_t(0)) + , dPdy_(measure_t(0)) + {} + + // image id: + unsigned int id_; + + // image value: + measure_t P_; + + // partial derivative with respect to x: + measure_t dPdx_; + + // partial derivative with respect to y: + measure_t dPdy_; + + // a pointer to the Jacobian of the transform: + const jacobian_t * J_; + }; + + // calculate the bounding box for a given set of image bounding boxes: + void + CalcMosaicBBox(point_t & mosaic_min, + point_t & mosaic_max, + std::vector & image_min, + std::vector & image_max) const; + + // calculate the mosaic variance image as well as maximum and mean variance: + typename image_t::Pointer + variance(measure_t & max_var, measure_t & avg_var) const; + + typename image_t::Pointer + variance(const typename TImage::SpacingType & sp, + const pnt2d_t & mosaic_min, + const pnt2d_t & mosaic_max, + measure_t & max_var, + measure_t & avg_var) const; + + typename image_t::Pointer + variance(const typename TImage::SpacingType & sp, + const pnt2d_t & mosaic_min, + const typename TImage::SizeType & sz, + measure_t & max_var, + measure_t & avg_var) const; + + typename image_t::Pointer + variance() const + { + measure_t max_var; + measure_t avg_var; + return variance(max_var, avg_var); + } + + // mosaic images: + std::vector image_; + std::vector mask_; + + // image transforms (cascaded): + mutable std::vector transform_; + + // image interpolators: + std::vector interpolator_; + + // image gradients: + std::vector gradient_; + + // adresses of the individual transform parameter indices within the + // concatenated parameter vector: + std::vector> address_; + + // number of shared/unique parameters per transform: + unsigned int n_shared_; + unsigned int n_unique_; + + // this mask defines which of the individual transform parameters are + // active and will be included into the metric parameter vector: + std::vector param_active_; + +protected: + ImageMosaicVarianceMetric() + : n_shared_(0) + , n_unique_(0) + , param_active_(0) + , m_NumberOfPixelsCounted(0) + {} + + virtual ~ImageMosaicVarianceMetric() {} + + // standard ITK debugging helper: + void + PrintSelf(std::ostream & os, Indent indent) const; + + // number of pixels in the overlapping regions of the mosaic: + mutable unsigned long m_NumberOfPixelsCounted; + +private: + // unimplemented on purpose: + ImageMosaicVarianceMetric(const Self &); + void + operator=(const Self &); + + typename image_t::RegionType m_MosaicRegion; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkImageMosaicVarianceMetric.hxx" +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __itkImageMosaicVarianceMetric_h diff --git a/include/itkImageMosaicVarianceMetric.hxx b/include/itkImageMosaicVarianceMetric.hxx new file mode 100644 index 0000000..f107910 --- /dev/null +++ b/include/itkImageMosaicVarianceMetric.hxx @@ -0,0 +1,751 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkImageMosaicVarianceMetric.txx +// Author : Pavel A. Koshevoy +// Created : 2005/06/22 17:19 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A metric class measuring mean pixel variance within +// the mosaic tile overlap regions. The metric derivative + +#ifndef _itkImageMosaicVarianceMetric_hxx +#define _itkImageMosaicVarianceMetric_hxx + +// local includes: +#include "itkIRCommon.h" +#include "itkIRUtils.h" +#include "IRTerminator.h" + +// system includes: +#include + +// ITK includes: +#include +#include + + +namespace itk +{ + +//---------------------------------------------------------------- +// setup_param_map +// +template +void +ImageMosaicVarianceMetric::setup_param_map(const std::vector & param_shared, + const std::vector & param_active) +{ + const unsigned int num_transforms = transform_.size(); + assert(num_transforms > 0); + + const unsigned int n_params = param_shared.size(); + assert(n_params > 0 && n_params == param_active.size() && n_params == transform_[0]->GetNumberOfParameters()); + + // save the active parameters mask: + param_active_ = param_active; + + // count the number of shared/unique params: + n_shared_ = 0; + n_unique_ = 0; + for (unsigned int i = 0; i < n_params; i++) + { + n_shared_ += param_active[i] && param_shared[i]; + n_unique_ += param_active[i] && !param_shared[i]; + } + + // build a mapping into the concatenated parameters vector, + // shared parameters will be stored at the head of the vector + // (to help with debugging) followed by the unique parameters: + resize(address_, num_transforms, n_params); + + unsigned int u_off = n_shared_; + for (unsigned int i = 0; i < num_transforms; i++) + { + unsigned int s_off = 0; + for (unsigned int j = 0; j < n_params; j++) + { + if (!param_active[j]) + continue; + + if (param_shared[j]) + { + address_[i][j] = s_off; + s_off++; + } + else + { + address_[i][j] = u_off; + u_off++; + } + + // FIXME: + // cout << "address[" << i << "][" << j << "] = " << address_[i][j] + // << endl; + } + } +} + +//---------------------------------------------------------------- +// GetTransformParameters +// +template +typename ImageMosaicVarianceMetric::params_t +ImageMosaicVarianceMetric::GetTransformParameters() const +{ + const unsigned int num_transforms = transform_.size(); + assert(num_transforms > 0); + + const unsigned int n_params = param_active_.size(); + assert(n_params > 0); + + params_t parameters(GetNumberOfParameters()); + for (unsigned int i = 0; i < num_transforms; i++) + { + params_t params = transform_[i]->GetParameters(); + for (unsigned int k = 0; k < n_params; k++) + { + if (!param_active_[k]) + continue; + parameters[address_[i][k]] = params[k]; + } + } + + return parameters; +} + +//---------------------------------------------------------------- +// SetTransformParameters +// +template +void +ImageMosaicVarianceMetric::SetTransformParameters(const params_t & parameters) const +{ + const unsigned int num_transforms = transform_.size(); + assert(num_transforms > 0); + + const unsigned int n_params = param_active_.size(); + assert(n_params > 0); + + // extract individual transform parameter values: + for (unsigned int i = 0; i < num_transforms; i++) + { + params_t params = transform_[i]->GetParameters(); + for (unsigned int k = 0; k < n_params; k++) + { + if (!param_active_[k]) + continue; + params[k] = parameters[address_[i][k]]; + } + transform_[i]->SetParameters(params); + } +} + +//---------------------------------------------------------------- +// GetNumberOfParameters +// +template +unsigned int +ImageMosaicVarianceMetric::GetNumberOfParameters() const +{ + const unsigned int num_transforms = transform_.size(); + return n_shared_ + n_unique_ * num_transforms; +} + +//---------------------------------------------------------------- +// Initialize +// +template +void +ImageMosaicVarianceMetric::Initialize() throw(ExceptionObject) +{ +#if 0 + cout << "sizeof(pixel_t) = " << sizeof(pixel_t) << endl + << "sizeof(scalar_t) = " << sizeof(scalar_t) << endl + << "sizeof(measure_t) = " << sizeof(measure_t) << endl + << "sizeof(gradient_pixel_t) = " << sizeof(gradient_pixel_t) << endl + << endl; +#endif + + const unsigned int num_transforms = address_.size(); + + if (num_transforms == 0) + { + itkExceptionMacro(<< "call setup_param_map first"); + } + + for (unsigned int i = 0; i < num_transforms; i++) + { + if (!transform_[i]) + { + itkExceptionMacro(<< "One of the transforms is missing"); + } + + if (!image_[i]) + { + itkExceptionMacro(<< "One of the images is missing"); + } + else if (image_[i]->GetSource()) + { + // if the image is provided by a source, update the source: + image_[i]->GetSource()->Update(); + } + } + + // setup the interpolators, calculate the gradient images: + interpolator_.resize(num_transforms); + gradient_.resize(num_transforms); + for (unsigned int i = 0; i < num_transforms; i++) + { + interpolator_[i] = interpolator_t::New(); + interpolator_[i]->SetInputImage(image_[i]); + + // calculate the image gradient: + typename gradient_filter_t::Pointer gradient_filter = gradient_filter_t::New(); + gradient_filter->SetInput(image_[i]); + + const typename image_t::SpacingType & spacing = image_[i]->GetSpacing(); + double maximum_spacing = 0.0; + for (unsigned int j = 0; j < ImageDimension; j++) + { + maximum_spacing = std::max(maximum_spacing, spacing[j]); + } + + gradient_filter->SetSigma(maximum_spacing); + gradient_filter->SetNormalizeAcrossScale(true); + + try + { + gradient_filter->Update(); + } + catch (itk::ExceptionObject & exception) + { + // oops: + cerr << "gradient filter threw an exception:" << endl << exception << endl; + throw exception; + } + + gradient_[i] = gradient_filter->GetOutput(); + } + + // If there are any observers on the metric, call them to give the + // user code a chance to set parameters on the metric: + this->InvokeEvent(InitializeEvent()); +} + +//---------------------------------------------------------------- +// GetValueAndDerivative +// +// FIXME: the following code assumes a 2D mosaic: +// +template +void +ImageMosaicVarianceMetric::GetValueAndDerivative(const params_t & parameters, + measure_t & measure, + derivative_t & derivative) const +{ + WRAP(itk_terminator_t terminator("ImageMosaicVarianceMetric::GetValueAndDerivative")); + + typedef typename image_t::IndexType index_t; + typedef typename image_t::SpacingType spacing_t; + typedef typename image_t::RegionType::SizeType imagesz_t; + typedef typename image_t::RegionType region_t; + + itkDebugMacro("GetValueAndDerivative( " << parameters << " ) "); + + // shortcuts: + const unsigned int num_images = interpolator_.size(); + if (num_images == 0) + { + itkExceptionMacro(<< "You forgot to call Initialize()"); + } + + // number of parameters per transform: + const unsigned int n_params = param_active_.size(); + if (n_params == 0) + { + itkExceptionMacro(<< "You forgot to call setup_param_map()"); + } + + // the derivative vector size: + const unsigned int n_concat = GetNumberOfParameters(); + + // compare the previous parameters to the next set of parameters: +#if 0 + { + params_t prev = GetTransformParameters(); + params_t diff = parameters; + cout << "0 diff: "; + for (unsigned int i = 0; i < n_concat; i++) + { + diff[i] -= prev[i]; + cout << setw(8 + 3) << diff[i]; + if (i == (n_concat - 1) / 2) cout << endl << "1 diff: "; + else if (i == (n_concat - 1)) cout << endl; + } + cout << endl; + } +#endif + +#if 0 + { + cout << "new parameters:\n" << parameters << endl; + for (unsigned int i = 0; i < transform_.size(); i++) + { + cout << i << ". " << transform_[i] << endl; + } + } +#endif + + // update the transforms: + SetTransformParameters(parameters); + + // calculate image bounding boxes in the mosaic space, as well as the + // mosaic bounding box: + pnt2d_t mosaic_min = pnt2d(std::numeric_limits::max(), std::numeric_limits::max()); + pnt2d_t mosaic_max = pnt2d(-mosaic_min[0], -mosaic_min[1]); + std::vector min(num_images); + std::vector max(num_images); + CalcMosaicBBox(mosaic_min, mosaic_max, min, max); + + // mosaic will have the same spacing as the first image in the mosaic: + spacing_t spacing = image_[0]->GetSpacing(); + imagesz_t mosaic_sz; + { + double nx = (mosaic_max[0] - mosaic_min[0]) / spacing[0]; + double ny = (mosaic_max[1] - mosaic_min[1]) / spacing[1]; + + mosaic_sz[0] = (unsigned int)(nx + 0.5); + mosaic_sz[1] = (unsigned int)(ny + 0.5); + } + + // update the region of interest if necessary: + region_t ROI = m_MosaicRegion; + if (ROI.GetNumberOfPixels() == 0) + { + // use largest possible region: + index_t index; + index[0] = 0; + index[1] = 0; + ROI.SetIndex(index); + ROI.SetSize(mosaic_sz); + } + + // setup the mosaic image, but don't allocate any buffers for it: + typename image_t::Pointer mosaic = image_t::New(); + mosaic->SetRegions(mosaic_sz); + mosaic->SetSpacing(spacing); + mosaic->SetOrigin(pnt2d(mosaic_min[0], mosaic_min[1])); + + // reset the accumulators: + derivative = derivative_t(n_concat); + derivative.Fill(NumericTraits::Zero); + measure = NumericTraits::Zero; + m_NumberOfPixelsCounted = 0; + + // pre-allocate and initialize image data for each image: + std::vector image_data(num_images); + for (unsigned int i = 0; i < num_images; i++) + { + image_data[i].id_ = i; + } + + // preallocate derivatives of the mean and variance: + std::vector dMu(n_concat); + std::vector dV(n_concat); + + // iterate over the mosaic, evaluate the metric in the overlapping regions: + typedef itk::ImageRegionConstIteratorWithIndex itex_t; + // bool FIXME_FIRST_RUN = true; + itex_t itex(mosaic, ROI); + for (itex.GoToBegin(); !itex.IsAtEnd(); ++itex) + { + // make sure there hasn't been an interrupt: + WRAP(terminator.terminate_on_request()); + + pnt2d_t point; + mosaic->TransformIndexToPhysicalPoint(itex.GetIndex(), point); + + // find pixels in overlapping regions of the mosaic: + std::list overlap; + for (unsigned int k = 0; k < num_images; k++) + { + if (point[0] < min[k][0] || point[0] > max[k][0] || point[1] < min[k][1] || point[1] > max[k][1]) + continue; + + typename transform_t::Pointer & t = transform_[k]; + pnt2d_t pt_k = t->TransformPoint(point); + index_t index; + // make sure the pixel maps into the image: + image_[k]->TransformPhysicalPointToIndex(pt_k, index); + if (!interpolator_[k]->IsInsideBuffer(pt_k)) + continue; + + // make sure the pixel maps into the mask: + if (mask_[k].GetPointer() && mask_[k]->GetPixel(index) == 0) + continue; + + // shortcut: + image_data_t & data = image_data[k]; + + // get the image value: + // data.P_ = image_[k]->GetPixel(index); // faster + data.P_ = interpolator_[k]->Evaluate(pt_k); // slower + + // get the image gradient: + const gradient_pixel_t & gradient = gradient_[k]->GetPixel(index); + data.dPdx_ = gradient[0]; + data.dPdy_ = gradient[1]; + + // get the transform Jacobian: + data.J_ = &(t->GetJacobian(point)); + + // add to the list: + overlap.push_back(&data); + } + + // skip over the regions that do not overlap: + if (overlap.size() < 2) + continue; + + // found a pixel in an overlapping region, increment the counter: + m_NumberOfPixelsCounted++; + + // shortcut: + measure_t normalization_factor = measure_t(1) / measure_t(overlap.size()); + + // calculate the mean and derivative of the mean: + measure_t Mu = measure_t(0); + dMu.assign(n_concat, measure_t(0)); + for (typename std::list::const_iterator i = overlap.begin(); i != overlap.end(); ++i) + { + const image_data_t & data = *(*i); + const unsigned int * addr = &(address_[data.id_][0]); + const jacobian_t & J = *(data.J_); + + Mu += data.P_; + + for (unsigned int j = 0; j < n_params; j++) + { + dMu[addr[j]] += (J[0][j] * data.dPdx_ + J[1][j] * data.dPdy_) * normalization_factor; + } + } + Mu *= normalization_factor; +#if 1 + // normalize the shared portion of dMu: + for (unsigned int i = 0; i < n_shared_; i++) + { + dMu[i] *= normalization_factor; + } +#endif + + // calculate the variance: + measure_t V = measure_t(0); + dV.assign(n_concat, measure_t(0)); + for (typename std::list::const_iterator i = overlap.begin(); i != overlap.end(); ++i) + { + const image_data_t & data = *(*i); + const unsigned int * addr = &(address_[data.id_][0]); + const jacobian_t & J = *(data.J_); + + measure_t d = data.P_ - Mu; + V += d * d; + + for (unsigned int j = 0; j < n_params; j++) + { + dV[addr[j]] += + measure_t(2) * d * (J[0][j] * data.dPdx_ + J[1][j] * data.dPdy_ - dMu[addr[j]]) * normalization_factor; + } + } + V *= normalization_factor; +#if 1 + // normalize the shared portion of dV: + for (unsigned int i = 0; i < n_shared_; i++) + { + dV[i] *= normalization_factor; + } +#endif + +#if 0 + if (FIXME_FIRST_RUN && overlap.size() > 2) + { + FIXME_FIRST_RUN = false; + for (typename std::list::const_iterator + i = overlap.begin(); i != overlap.end(); ++i) + { + const image_data_t & data = *(*i); + const unsigned int * addr = &(address_[data.id_][0]); + const jacobian_t & J = *(data.J_); + + for (unsigned int k = 0; k < 2; k++) + { + cout << "J[" << k << "] ="; + for (unsigned int j = 0; j < n_params; j++) + { + cout << ' ' << J[k][j]; + } + cout << endl; + } + cout << endl; + } + } +#endif + // update the measure: + measure += V; + + // update the derivative: + for (unsigned int i = 0; i < n_concat; i++) + { + derivative[i] += dV[i]; + } + } + + if (m_NumberOfPixelsCounted == 0) + { + itkExceptionMacro(<< "mosaic contains no overlapping images"); + return; + } + + measure_t normalization_factor = measure_t(1) / measure_t(m_NumberOfPixelsCounted); + measure *= normalization_factor; + + for (unsigned int i = 0; i < n_concat; i++) + { + derivative[i] *= normalization_factor; + } + +#if 0 + cout << "stdV: " << measure << endl; + cout << "0 dV: "; + for (unsigned int i = 0; i < n_concat; i++) + { + cout << setw(7 + 3) << derivative[i]; + if (i == (n_concat - 1) / 2) cout << endl << "1 dV: "; + else if (i == (n_concat - 1)) cout << endl; + else cout << ' '; + } + cout << endl; +#endif +} + +//---------------------------------------------------------------- +// CalcMosaicBBox +// +template +void +ImageMosaicVarianceMetric::CalcMosaicBBox(point_t & mosaic_min, + point_t & mosaic_max, + std::vector & min, + std::vector & max) const +{ + // image space bounding boxes: + std::vector image_min; + std::vector image_max; + calc_image_bboxes(image_, image_min, image_max); + + std::vector ConstTransforms; + ConstTransforms.assign(transform_.begin(), transform_.end()); + + // mosaic space bounding boxes: + calc_mosaic_bboxes(ConstTransforms, image_min, image_max, min, max); + + // mosiac bounding box: + calc_mosaic_bbox(min, max, mosaic_min, mosaic_max); +} + +//---------------------------------------------------------------- +// variance +// +template +typename TImage::Pointer +ImageMosaicVarianceMetric::variance(measure_t & max_var, measure_t & avg_var) const +{ + max_var = measure_t(0); + avg_var = measure_t(0); + + // shortcut: + const unsigned int num_images = interpolator_.size(); + if (num_images == 0) + { + itkExceptionMacro(<< "You forgot to call Initialize()"); + } + + pnt2d_t mosaic_min = pnt2d(std::numeric_limits::max(), std::numeric_limits::max()); + pnt2d_t mosaic_max = pnt2d(-mosaic_min[0], -mosaic_min[1]); + std::vector min(num_images); + std::vector max(num_images); + CalcMosaicBBox(mosaic_min, mosaic_max, min, max); + + return variance(image_[0]->GetSpacing(), mosaic_min, mosaic_max, max_var, avg_var); +} + +//---------------------------------------------------------------- +// variance +// +template +typename TImage::Pointer +ImageMosaicVarianceMetric::variance(const typename TImage::SpacingType & mosaic_sp, + const pnt2d_t & mosaic_min, + const pnt2d_t & mosaic_max, + measure_t & max_var, + measure_t & avg_var) const +{ + typename TImage::SizeType mosaic_sz; + mosaic_sz[0] = (unsigned int)((mosaic_max[0] - mosaic_min[0]) / mosaic_sp[0]); + mosaic_sz[1] = (unsigned int)((mosaic_max[1] - mosaic_min[1]) / mosaic_sp[1]); + return variance(mosaic_sp, mosaic_min, mosaic_sz, max_var, avg_var); +} + +//---------------------------------------------------------------- +// variance +// +template +typename TImage::Pointer +ImageMosaicVarianceMetric::variance(const typename TImage::SpacingType & mosaic_sp, + const pnt2d_t & mosaic_min, + const typename TImage::SizeType & mosaic_sz, + measure_t & max_var, + measure_t & avg_var) const +{ + WRAP(itk_terminator_t terminator("ImageMosaicVarianceMetric::GetValueAndDerivative")); + + typedef typename image_t::IndexType index_t; + typedef typename image_t::RegionType::SizeType imagesz_t; + + max_var = measure_t(0); + avg_var = measure_t(0); + + // shortcut: + const unsigned int num_images = interpolator_.size(); + if (num_images == 0) + { + itkExceptionMacro(<< "You forgot to call Initialize()"); + } + + // calculate image bounding boxes in the mosaic space, as well as the + // mosaic bounding box: + pnt2d_t global_min = pnt2d(std::numeric_limits::max(), std::numeric_limits::max()); + pnt2d_t global_max = pnt2d(-global_min[0], -global_min[1]); + std::vector min(num_images); + std::vector max(num_images); + CalcMosaicBBox(global_min, global_max, min, max); + + // setup the mosaic image: + typename image_t::Pointer mosaic = image_t::New(); + mosaic->SetOrigin(mosaic_min); + mosaic->SetRegions(mosaic_sz); + mosaic->SetSpacing(mosaic_sp); + mosaic->Allocate(); + mosaic->FillBuffer(NumericTraits::Zero); + + // iterate over the mosaic, evaluate the metric in the overlapping regions: + typedef itk::ImageRegionIteratorWithIndex itex_t; + unsigned int pixel_count = 0; + + itex_t itex(mosaic, mosaic->GetLargestPossibleRegion()); + for (itex.GoToBegin(); !itex.IsAtEnd(); ++itex) + { + // make sure there hasn't been an interrupt: + WRAP(terminator.terminate_on_request()); + + pnt2d_t point; + mosaic->TransformIndexToPhysicalPoint(itex.GetIndex(), point); + + // find pixels in overlapping regions of the mosaic: + std::list overlap; + for (unsigned int k = 0; k < num_images; k++) + { + if (point[0] < min[k][0] || point[0] > max[k][0] || point[1] < min[k][1] || point[1] > max[k][1]) + continue; + + const typename transform_t::Pointer & t = transform_[k]; + pnt2d_t pt_k = t->TransformPoint(point); + index_t index; + image_[k]->TransformPhysicalPointToIndex(pt_k, index); + + // make sure the pixel maps into the image: + if (!interpolator_[k]->IsInsideBuffer(pt_k)) + continue; + + // make sure the pixel maps into the mask: + if (mask_[k].GetPointer() && mask_[k]->GetPixel(index) == 0) + continue; + + // get the image value, add it to the list: + overlap.push_back(interpolator_[k]->Evaluate(pt_k)); + } + + // skip over the regions that do not overlap: + if (overlap.size() < 2) + continue; + + // found a pixel in an overlapping region, increment the counter: + pixel_count++; + + // shortcut: + measure_t normalization_factor = measure_t(1) / measure_t(overlap.size()); + + // calculate the mean and derivative of the mean: + measure_t Mu = measure_t(0); + for (std::list::const_iterator i = overlap.begin(); i != overlap.end(); ++i) + { + Mu += *i; + } + Mu *= normalization_factor; + + // calculate the variance: + measure_t V = measure_t(0); + for (std::list::const_iterator i = overlap.begin(); i != overlap.end(); ++i) + { + measure_t d = *i - Mu; + V += d * d; + } + V *= normalization_factor; + itex.Set(pixel_t(V)); + + max_var = std::max(max_var, V); + avg_var += V; + } + + measure_t normalization_factor = measure_t(1) / measure_t(pixel_count); + avg_var *= normalization_factor; + + return mosaic; +} + +//---------------------------------------------------------------- +// PrintSelf +// +template +void +ImageMosaicVarianceMetric::PrintSelf(std::ostream & os, Indent indent) const +{ + // FIXME: write me: + + Superclass::PrintSelf(os, indent); + os << indent << "MosaicRegion: " << m_MosaicRegion << std::endl + << indent << "Number of Pixels Counted: " << m_NumberOfPixelsCounted << std::endl; +} + +} // end namespace itk + + +#endif // _itkImageMosaicVarianceMetric_txx diff --git a/include/itkMeshTransform.h b/include/itkMeshTransform.h new file mode 100644 index 0000000..be1e65b --- /dev/null +++ b/include/itkMeshTransform.h @@ -0,0 +1,356 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: nil -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkMeshTransform.h +// Author : Pavel A. Koshevoy +// Created : 2009/05/06 20:13 +// Copyright : (C) 2009 University of Utah +// License : GPLv2 +// Description : A discontinuous transform -- a Delaunay triangle mesh +// mapped to an image. At each vertex, in addition to image +// space coordinates a second set of coordinates is stored. +// This is similar to texture mapped OpenGL triangle meshes, +// where the texture coordinates correspond to the image space +// vertex coordinates. + +#ifndef __itkMeshTransform_h +#define __itkMeshTransform_h + +// system includes: +#include +#include +#include + +// ITK includes: +#include +#include +#include +#include +#include + +// local includes: +#include "IRGridTransform.h" + + +//---------------------------------------------------------------- +// itk::MeshTransform +// +namespace itk +{ +class MeshTransform : public Transform +{ +public: + // standard typedefs: + typedef MeshTransform Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + typedef Transform Superclass; + + // Base inverse transform type: + typedef Superclass::InverseTransformType InverseTransformType; + typedef SmartPointer InverseTransformPointer; + + // RTTI: + itkTypeMacro(MeshTransform, Transform); + + // macro for instantiation through the object factory: + itkNewMacro(Self); + + /** Standard scalar type for this class. */ + typedef double ScalarType; + + /** Dimension of the domain space. */ + itkStaticConstMacro(InputSpaceDimension, unsigned int, 2); + itkStaticConstMacro(OutputSpaceDimension, unsigned int, 2); + + // shortcuts: + typedef Superclass::ParametersType ParametersType; + typedef Superclass::JacobianType JacobianType; + + typedef Superclass::InputPointType InputPointType; + typedef Superclass::OutputPointType OutputPointType; + + // virtual: + OutputPointType + TransformPoint(const InputPointType & x) const + { + OutputPointType y; + if (transform_.grid_.cols_ == 0 || transform_.grid_.cols_ == 0) + { + // No grid has been setup. Use the identity. + y[0] = x[0]; + y[1] = x[1]; + } + else if (is_inverse()) + { + pnt2d_t uv; + uv[0] = (x[0] - transform_.tile_min_[0]) / transform_.tile_ext_[0]; + uv[1] = (x[1] - transform_.tile_min_[1]) / transform_.tile_ext_[1]; + transform_.transform_inv(uv, y); + } + else + { + transform_.transform(x, y); + y[0] *= transform_.tile_ext_[0]; + y[1] *= transform_.tile_ext_[1]; + y[0] += transform_.tile_min_[0]; + y[1] += transform_.tile_min_[1]; + } + + // ITK is not forgiving to NaNs: + if (y[0] != y[0]) + { + y[0] = std::numeric_limits::max(); + y[1] = y[0]; + } + + return y; + } + + // Inverse transformations: + // If y = Transform(x), then x = BackTransform(y); + // if no mapping from y to x exists, then an exception is thrown. + InputPointType + BackTransformPoint(const OutputPointType & y) const + { + InputPointType x; + if (is_inverse()) + { + transform_.transform(y, x); + x[0] *= transform_.tile_ext_[0]; + x[1] *= transform_.tile_ext_[1]; + x[0] += transform_.tile_min_[0]; + x[1] += transform_.tile_min_[1]; + } + else + { + pnt2d_t uv; + uv[0] = (y[0] - transform_.tile_min_[0]) / transform_.tile_ext_[0]; + uv[1] = (y[1] - transform_.tile_min_[1]) / transform_.tile_ext_[1]; + transform_.transform_inv(uv, x); + } + + // ITK does not handle NaNs well: + if (x[0] != x[0]) + { + x[0] = std::numeric_limits::max(); + x[1] = x[0]; + } + + return x; + } + + // virtual: + void + SetFixedParameters(const ParametersType & params) + { + this->m_FixedParameters = params; + } + + // virtual: + const ParametersType & + GetFixedParameters() const + { + ParametersType params = this->m_FixedParameters; + + // acceleration grid size: + params[1] = transform_.grid_.rows_; + params[2] = transform_.grid_.cols_; + + // bounding box if the image associated with this transform: + params[3] = transform_.tile_min_[0]; + params[4] = transform_.tile_min_[1]; + params[5] = transform_.tile_ext_[0]; + params[6] = transform_.tile_ext_[1]; + + // number of vertices in the Delaunay triangle mesh: + params[7] = transform_.grid_.mesh_.size(); + + // update the parameters vector: + MeshTransform * fake = const_cast(this); + fake->m_FixedParameters = params; + return this->m_FixedParameters; + } + + // virtual: + void + SetParameters(const ParametersType & params) + { + this->m_Parameters = params; + + std::size_t accel_grid_rows = std::size_t(this->m_FixedParameters[1]); + std::size_t accel_grid_cols = std::size_t(this->m_FixedParameters[2]); + + pnt2d_t tile_min; + tile_min[0] = this->m_FixedParameters[3]; + tile_min[1] = this->m_FixedParameters[4]; + + pnt2d_t tile_max; + tile_max[0] = tile_min[0] + this->m_FixedParameters[5]; + tile_max[1] = tile_min[1] + this->m_FixedParameters[6]; + + const std::size_t num_points = std::size_t(this->m_FixedParameters[7]); + std::vector uv(num_points); + std::vector xy(num_points); + + for (unsigned int i = 0; i < num_points; i++) + { + const std::size_t idx = i * 4; + + uv[i][0] = params[idx + 0]; + uv[i][1] = params[idx + 1]; + xy[i][0] = params[idx + 2]; + xy[i][1] = params[idx + 3]; + } + + transform_.setup(tile_min, tile_max, uv, xy, accel_grid_rows, accel_grid_cols); + } + + // virtual: + const ParametersType & + GetParameters() const + { + ParametersType params(GetNumberOfParameters()); + + const std::size_t num_points = transform_.grid_.mesh_.size(); + for (std::size_t i = 0; i < num_points; i++) + { + const vertex_t & vx = transform_.grid_.mesh_[i]; + const std::size_t idx = i * 4; + + params[idx + 0] = vx.uv_[0]; + params[idx + 1] = vx.uv_[1]; + params[idx + 2] = vx.xy_[0]; + params[idx + 3] = vx.xy_[1]; + } + + // update the parameters vector: + MeshTransform * fake = const_cast(this); + fake->m_Parameters = params; + return this->m_Parameters; + } + + // virtual: + unsigned int + GetNumberOfParameters() const + { + return 4 * transform_.grid_.mesh_.size(); + } + + // virtual: return an inverse of this transform. + InverseTransformPointer + GetInverse() const + { + MeshTransform::Pointer inv = MeshTransform::New(); + inv->setup(transform_, !is_inverse()); + return inv.GetPointer(); + } + + // setup the transform: + void + setup(const the_mesh_transform_t & transform, const bool & is_inverse = false) + { + transform_ = transform; + GetParameters(); + GetFixedParameters(); + this->m_Jacobian.SetSize(2, GetNumberOfParameters()); + this->m_FixedParameters[0] = is_inverse ? 1.0 : 0.0; + } + + // inverse transform flag check: + inline bool + is_inverse() const + { + return this->m_FixedParameters[0] != 0.0; + } + + // virtual: + const JacobianType & + GetJacobian(const InputPointType & point) const + { + // FIXME: 20061227 -- this function was written and not tested: + + // these scales are necessary to account for the fact that + // the_mesh_transform_t expects uv in the [0,1]x[0,1] range, + // where as we remap it into the image tile physical coordinates + // according to tile_min_ and tile_ext_: + double Su = transform_.tile_ext_[0]; + double Sv = transform_.tile_ext_[1]; + + unsigned int idx[3]; + double jac[12]; + this->m_Jacobian.SetSize(2, GetNumberOfParameters()); + this->m_Jacobian.Fill(0.0); + if (transform_.jacobian(point, idx, jac)) + { + for (unsigned int i = 0; i < 3; i++) + { + unsigned int addr = idx[i] * 2; + this->m_Jacobian(0, addr) = Su * jac[i * 2]; + this->m_Jacobian(0, addr + 1) = Su * jac[i * 2 + 1]; + this->m_Jacobian(1, addr) = Sv * jac[i * 2 + 6]; + this->m_Jacobian(1, addr + 1) = Sv * jac[i * 2 + 7]; + } + } + + return this->m_Jacobian; + } + +protected: + MeshTransform() + : Superclass(2, 0) + { + this->m_FixedParameters.SetSize(8); + + // initialize the inverse flag: + this->m_FixedParameters[0] = 0.0; + + // acceleration grid size: + this->m_FixedParameters[1] = 0.0; + this->m_FixedParameters[2] = 0.0; + + // bounding box if the image associated with this transform: + this->m_FixedParameters[3] = std::numeric_limits::max(); + this->m_FixedParameters[4] = this->m_FixedParameters[3]; + this->m_FixedParameters[5] = -(this->m_FixedParameters[3]); + this->m_FixedParameters[6] = -(this->m_FixedParameters[3]); + + // number of vertices in the Delaunay triangle mesh: + this->m_FixedParameters[7] = 0.0; + } + +private: + // disable default copy constructor and assignment operator: + MeshTransform(const Self & other); + const Self & + operator=(const Self & t); + +public: + // the actual transform: + the_mesh_transform_t transform_; + +}; // class MeshTransform + +} // namespace itk + + +#endif // __itkMeshTransform_h diff --git a/include/itkRadialDistortionTransform.h b/include/itkRadialDistortionTransform.h new file mode 100644 index 0000000..d2b5e90 --- /dev/null +++ b/include/itkRadialDistortionTransform.h @@ -0,0 +1,255 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkRadialDistortionTransform.h +// Author : Pavel A. Koshevoy +// Created : 2005/06/03 10:16 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A radial distortion transform. + +#ifndef __itkRadialDistortionTransform_h +#define __itkRadialDistortionTransform_h + +// system includes: +#include +#include + +// ITK includes: +#include +#include +#include + +// local includes: +#include "itkInverseTransform.h" + + +//---------------------------------------------------------------- +// itk::RadialDistortionTransform +// +// Let +// A = (a + ta * Rmax - ac) +// B = (b + tb * Rmax - bc) +// where ta, tb specify percentage of translation along the +// a and b axis respectively, and +// ac = 0.5 * (a_min + a_max) +// bc = 0.5 * (b_min + b_max) +// define the center of distortion (specified by a_min, a_max, +// b_min, b_max -- the bounding box of some image). +// +// The transform is defined as +// x(a, b) = ac + A * S +// y(a, b) = bc + B * S +// +// where +// S = sum(n in [0, N - 1], k[n] * pow(R2 / Rmax2, n)); +// R2 = A^2 + B^2 +// Rmax2 = Rmax * Rmax +// +namespace itk +{ +template +class RadialDistortionTransform : public Transform +{ +public: + // standard typedefs: + typedef RadialDistortionTransform Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + typedef Transform Superclass; + + // Base inverse transform type: + typedef typename Superclass::InverseTransformType InverseTransformType; + typedef SmartPointer InverseTransformPointer; + + // static constant for the number of polynomial coefficients: + itkStaticConstMacro(Nk, unsigned int, N); + + // RTTI: + itkTypeMacro(RadialDistortionTransform, Transform); + + // macro for instantiation through the object factory: + itkNewMacro(Self); + + /** Standard scalar type for this class. */ + typedef typename Superclass::ScalarType ScalarType; + + // shortcuts: + typedef typename Superclass::ParametersType ParametersType; + typedef typename Superclass::JacobianType JacobianType; + + typedef typename Superclass::InputPointType InputPointType; + typedef typename Superclass::OutputPointType OutputPointType; + + // virtual: + OutputPointType + TransformPoint(const InputPointType & p) const; + + // virtual: Inverse transformations: + // If y = Transform(x), then x = BackTransform(y); + // if no mapping from y to x exists, then an exception is thrown. + InputPointType + BackTransformPoint(const OutputPointType & y) const; + + // virtual: + void + SetFixedParameters(const ParametersType & params) + { + this->m_FixedParameters = params; + } + + // virtual: + const ParametersType & + GetFixedParameters() const + { + return this->m_FixedParameters; + } + + // virtual: + void + SetParameters(const ParametersType & params) + { + this->m_Parameters = params; + } + + // virtual: + const ParametersType & + GetParameters() const + { + return this->m_Parameters; + } + + // virtual: mumber of parameters that define this transform: + unsigned int + GetNumberOfParameters() const + { + return N + 2; + } + + // virtual: + const JacobianType & + GetJacobian(const InputPointType & point) const; + + // virtual: return an inverse of this transform. + InverseTransformPointer + GetInverse() const + { + typedef InverseTransform InvTransformType; + typename InvTransformType::Pointer inv = InvTransformType::New(); + inv->SetForwardTransform(this); + return inv.GetPointer(); + } + + // setup the fixed transform parameters: + void + setup( // image bounding box expressed in the physical space: + const double a_min, + const double a_max, + const double b_min, + const double b_max, + + // normalization parameter (image radius in physical space): + const double Rmax = 0.0) + { + double & ac_ = this->m_FixedParameters[0]; + double & bc_ = this->m_FixedParameters[1]; + double & rmax_ = this->m_FixedParameters[2]; + + // center of distortion: + ac_ = 0.5 * (a_min + a_max); + bc_ = 0.5 * (b_min + b_max); + + // setup the normalization parameter: + if (Rmax != 0.0) + { + rmax_ = Rmax; + } + else + { + const double w = a_max - a_min; + const double h = b_max - b_min; + rmax_ = sqrt(w * w + h * h) / 2.0; + } + } + + // setup the translation parameters: + void + setup_translation( // translation is expressed in the physical space: + const double ta_Rmax = 0.0, + const double tb_Rmax = 0.0) + { + const double & Rmax = this->m_FixedParameters[2]; + assert(Rmax != 0.0); + + // store ta,tb as translation relative to Rmax: + double & ta = this->m_Parameters[N]; + double & tb = this->m_Parameters[N + 1]; + ta = ta_Rmax / Rmax; + tb = tb_Rmax / Rmax; + } + + // helper required by BackTransform: + // evaluate F = T(x), J = dT/dx (another Jacobian): + void + eval(const std::vector & x, std::vector & F, std::vector> & J) const; + + // accessors to the fixed normalization parameter: + inline const double & + GetRmax() const + { + return this->m_FixedParameters[2]; + } + + // generate a mask of shared parameters: + static void + setup_shared_params_mask(bool shared, std::vector & mask) + { + mask.assign(N + 2, false); + for (unsigned int i = 0; i < N; i++) + { + mask[i] = shared; + } + } + +protected: + RadialDistortionTransform(); + + // virtual: + void + PrintSelf(std::ostream & s, Indent indent) const; + +private: + // disable default copy constructor and assignment operator: + RadialDistortionTransform(const Self & other); + const Self & + operator=(const Self & t); + +}; // class RadialDistortionTransform + +} // namespace itk + + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkRadialDistortionTransform.hxx" +#endif + +#endif // __itkRadialDistortionTransform_h diff --git a/include/itkRadialDistortionTransform.hxx b/include/itkRadialDistortionTransform.hxx new file mode 100644 index 0000000..5784f7c --- /dev/null +++ b/include/itkRadialDistortionTransform.hxx @@ -0,0 +1,282 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkRadialDistortionTransform.txx +// Author : Pavel A. Koshevoy +// Created : 2005/06/03 10:16 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A radial distortion transform. + +#ifndef _itkRadialDistortionTransform_txx +#define _itkRadialDistortionTransform_txx + +// local includes: +#include "itkNumericInverse.h" + +// system includes: +#include +#include + + +namespace itk +{ +//---------------------------------------------------------------- +// RadialDistortionTransform +// +template +RadialDistortionTransform::RadialDistortionTransform() + : Superclass(2, N + 2) // there are k0, k1,.. kN-1, tx, ty parameters: +{ + // by default, the parameters are initialized for an identity transform: + this->m_Parameters[0] = 1.0; + for (unsigned int i = 1; i < N + 2; i++) + { + this->m_Parameters[i] = 0.0; + } + + // allocate some space for the fixed parameters (ac, bc, Rmax): + this->m_FixedParameters.SetSize(3); + double & ac = this->m_FixedParameters[0]; + double & bc = this->m_FixedParameters[1]; + double & Rmax = this->m_FixedParameters[2]; + ac = 0.0; + bc = 0.0; + Rmax = 0.0; +} + +//---------------------------------------------------------------- +// TransformPoint +// +template +typename RadialDistortionTransform::OutputPointType +RadialDistortionTransform::TransformPoint(const InputPointType & x) const +{ + const double & ac = this->m_FixedParameters[0]; + const double & bc = this->m_FixedParameters[1]; + const double & Rmax = this->m_FixedParameters[2]; + + const double & ta = this->m_Parameters[N]; + const double & tb = this->m_Parameters[N + 1]; + + const ScalarType & a = x[0]; + const ScalarType & b = x[1]; + + const double A = (a + ta * Rmax - ac); + const double B = (b + tb * Rmax - bc); + const double A2 = A * A; + const double B2 = B * B; + const double R2 = A2 + B2; + const double Rmax2 = Rmax * Rmax; + const double RRmax2 = R2 / Rmax2; + + // n = 0: + double S = this->m_Parameters[0]; + + // n = 1 ... N-1: + double RRmax2n = RRmax2; + for (unsigned int n = 1; n < N; n++) + { + const double knRRmax2n = this->m_Parameters[n] * RRmax2n; + S += knRRmax2n; + RRmax2n *= RRmax2; + } + + OutputPointType y; + y[0] = ac + A * S; + y[1] = bc + B * S; + + return y; +} + +//---------------------------------------------------------------- +// BackTransformPoint +// +template +typename RadialDistortionTransform::InputPointType +RadialDistortionTransform::BackTransformPoint(const OutputPointType & y) const +{ + NumericInverse> inverse(*this); + + std::vector vy(2); + std::vector vx(2); + vy[0] = y[0]; + vy[1] = y[1]; + + // initialize x: first guess - x is close to y: + vx = vy; + const double & Rmax = this->m_FixedParameters[2]; + const double & ta = this->m_Parameters[N]; + const double & tb = this->m_Parameters[N + 1]; + vx[0] -= ta * Rmax; + vx[1] -= tb * Rmax; + bool ok = inverse.transform(vy, vx, true); + if (!ok) + assert(false); + + OutputPointType x; + x[0] = vx[0]; + x[1] = vx[1]; + + return x; +} + +//---------------------------------------------------------------- +// GetJacobian +// +template +const typename RadialDistortionTransform::JacobianType & +RadialDistortionTransform::GetJacobian(const InputPointType & x) const +{ + const double & ac = this->m_FixedParameters[0]; + const double & bc = this->m_FixedParameters[1]; + const double & Rmax = this->m_FixedParameters[2]; + + const double & ta = this->m_Parameters[N]; + const double & tb = this->m_Parameters[N + 1]; + + const ScalarType & a = x[0]; + const ScalarType & b = x[1]; + + const double A = (a + ta * Rmax - ac); + const double B = (b + tb * Rmax - bc); + const double A2 = A * A; + const double B2 = B * B; + const double R2 = A2 + B2; + const double Rmax2 = Rmax * Rmax; + const double RRmax2 = R2 / Rmax2; + + // derivatives with respect to kn: + double RRmax2n = double(1); + for (unsigned int n = 0; n < N; n++) + { + this->m_Jacobian(0, n) = A * RRmax2n; + this->m_Jacobian(1, n) = B * RRmax2n; + RRmax2n *= RRmax2; + } + + // derivatives with respect to ta, tb: + + // n = 0: + double P = double(0); + double Q = double(0); + + // n = 1 ... N-1: + double RRmax2n1 = double(1); // (R^2 / Rmax^2)^(n - 1) + for (unsigned int n = 1; n < N; n++) + { + const double & k = this->m_Parameters[n]; + const double scale_n = k * RRmax2n1; // kn * (R^2/Rmax^2)^n + P += scale_n; + Q += scale_n * double(n); + RRmax2n1 *= RRmax2; + } + + double S = this->m_Parameters[0] + RRmax2 * P; + + this->m_Jacobian(0, N) = Rmax * S + (double(2) * A2 / Rmax) * Q; // dx/dta + this->m_Jacobian(1, N + 1) = Rmax * S + (double(2) * B2 / Rmax) * Q; // dy/dtb + + this->m_Jacobian(0, N + 1) = ((double(2) * A * B) / Rmax) * Q; // dx/dtb + this->m_Jacobian(1, N) = this->m_Jacobian(0, N + 1); // dy/dta + + return this->m_Jacobian; +} + +//---------------------------------------------------------------- +// eval +// +template +void +RadialDistortionTransform::eval(const std::vector & x, + std::vector & F, + std::vector> & J) const +{ + const double & ac = this->m_FixedParameters[0]; + const double & bc = this->m_FixedParameters[1]; + const double & Rmax = this->m_FixedParameters[2]; + + const double & ta = this->m_Parameters[N]; + const double & tb = this->m_Parameters[N + 1]; + + const ScalarType & a = x[0]; + const ScalarType & b = x[1]; + + const double A = (a + ta * Rmax - ac); + const double B = (b + tb * Rmax - bc); + const double A2 = A * A; + const double B2 = B * B; + const double R2 = A2 + B2; + const double Rmax2 = Rmax * Rmax; + const double RRmax2 = R2 / Rmax2; + + // n = 0: + double P = double(0); + double Q = double(0); + + // n = 1 ... N-1: + double RRmax2n1 = double(1); // (R^2 / Rmax^2)^(n - 1) + for (unsigned int n = 1; n < N; n++) + { + const double & k = this->m_Parameters[n]; + const double scale_n = k * RRmax2n1; // kn * (R^2/Rmax^2)^n + P += scale_n; + Q += scale_n * double(n); + RRmax2n1 *= RRmax2; + } + + double S = this->m_Parameters[0] + RRmax2 * P; + F[0] = ac + A * S; + F[1] = bc + B * S; + + // calc dF/dx: + J[0][0] = Rmax * S + (double(2) * A2 / Rmax) * Q; + J[1][1] = Rmax * S + (double(2) * B2 / Rmax) * Q; + + J[0][1] = ((double(2) * A * B) / Rmax) * Q; + J[1][0] = J[0][1]; +} + +//---------------------------------------------------------------- +// PrintSelf +// +template +void +RadialDistortionTransform::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + for (unsigned int i = 0; i < N; i++) + { + os << indent << "k[" << i << "] = " << this->m_Parameters[i] << std::endl; + } + + os << indent << "ta = " << this->m_Parameters[N] << std::endl + << indent << "tb = " << this->m_Parameters[N + 1] << std::endl + << indent << "ac = " << this->m_FixedParameters[0] << std::endl + << indent << "bc = " << this->m_FixedParameters[1] << std::endl + << indent << "Rmax = " << this->m_FixedParameters[2] << std::endl; +} + +} // namespace itk + + +#endif // _itkRadialDistortionTransform_txx diff --git a/include/itkRegularStepGradientDescentOptimizer2.h b/include/itkRegularStepGradientDescentOptimizer2.h new file mode 100644 index 0000000..44b6ce6 --- /dev/null +++ b/include/itkRegularStepGradientDescentOptimizer2.h @@ -0,0 +1,196 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkRegularStepGradientDescentOptimizer2.h +// Author : Pavel A. Koshevoy, Tolga Tasdizen +// Created : 2005/11/11 14:54 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An enhanced version of the +// itk::RegularStepGradientDescentOptimizer +// fixing a bug with relaxation and adding support for +// step size increase (pick up pace), back tracking and +// keeping track of the best metric value +// and associated parameters. + +#ifndef __itkRegularStepGradientDescentOptimizer2_h +#define __itkRegularStepGradientDescentOptimizer2_h + +// ITK includes: +#include + +// the includes: +#include + + +namespace itk +{ +#if 0 +} +#endif + +//---------------------------------------------------------------- +// RegularStepGradientDescentOptimizer2 +// +class ITK_EXPORT RegularStepGradientDescentOptimizer2 : + public SingleValuedNonLinearOptimizer +{ +public: + /** Standard "Self" typedef. */ + typedef RegularStepGradientDescentOptimizer2 Self; + typedef SingleValuedNonLinearOptimizer Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(RegularStepGradientDescentOptimizer2, + SingleValuedNonLinearOptimizer); + + //---------------------------------------------------------------- + // StopConditionType + // + // Codes of stopping conditions: + // + typedef enum { + GradientMagnitudeTolerance = 1, + StepTooSmall = 2, + ImageNotAvailable = 3, + SamplesNotAvailable = 4, + MaximumNumberOfIterations = 5 + } StopConditionType; + + /** Specify whether to minimize or maximize the cost function. */ + itkSetMacro(Maximize, bool); + itkGetConstReferenceMacro(Maximize, bool); + itkBooleanMacro(Maximize); + + inline bool GetMinimize() const { return !m_Maximize; } + inline void SetMinimize(bool v) { this->SetMaximize(!v); } + inline void MinimizeOn() { SetMaximize(false); } + inline void MinimizeOff() { SetMaximize(true); } + + // virtual: + void StartOptimization(); + + // virtual: + void ResumeOptimization(); + + // virtual: + void StopOptimization(); + + // Set/Get parameters to control the optimization process: + itkSetMacro(MaximumStepLength, double); + itkSetMacro(MinimumStepLength, double); + itkGetConstReferenceMacro(MaximumStepLength, double); + itkGetConstReferenceMacro(MinimumStepLength, double); + + itkSetMacro(RelaxationFactor, double); + itkGetConstReferenceMacro(RelaxationFactor, double); + + itkSetMacro(NumberOfIterations, unsigned long); + itkGetConstReferenceMacro(NumberOfIterations, unsigned long); + + itkSetMacro(GradientMagnitudeTolerance, double); + itkGetConstReferenceMacro(GradientMagnitudeTolerance, double); + + itkSetMacro(BackTracking, bool); + itkGetConstReferenceMacro(BackTracking, bool); + + itkSetMacro(PickUpPaceSteps, unsigned int); + itkGetConstReferenceMacro(PickUpPaceSteps, unsigned int); + + inline void SetLog(the_log_t * log) + { log_ = log; } + + itkGetConstReferenceMacro(CurrentStepLength, double); + itkGetConstMacro(CurrentIteration, unsigned int); + itkGetConstReferenceMacro(StopCondition, StopConditionType); + itkGetConstReferenceMacro(Value, MeasureType); + itkGetConstReferenceMacro(Gradient, DerivativeType); + itkGetConstReferenceMacro(BestParams, ParametersType); + itkGetConstReferenceMacro(BestValue, MeasureType); + +protected: + RegularStepGradientDescentOptimizer2(); + virtual ~RegularStepGradientDescentOptimizer2() {} + + // advance one step following the gradient direction: + virtual void AdvanceOneStep(); + + // advance one step along the given direction vector + // scaled by the step length scale; this method is + // called by AdvanceOneStep: + virtual void StepAlongGradient(double step_length_scale, + const DerivativeType & step_direction); + + // virtual: + void PrintSelf(std::ostream& os, Indent indent) const; + + // data members: + DerivativeType m_Gradient; + DerivativeType m_PreviousGradient; + + bool m_Stop; + bool m_Maximize; + MeasureType m_Value; + MeasureType m_PreviousValue; + double m_GradientMagnitudeTolerance; + double m_MaximumStepLength; + double m_MinimumStepLength; + double m_CurrentStepLength; + double m_RelaxationFactor; + StopConditionType m_StopCondition; + unsigned long m_NumberOfIterations; + unsigned long m_CurrentIteration; + + // this keeps track of the best function value and corresponding parameters: + MeasureType m_BestValue; + ParametersType m_BestParams; + + // this flag controls whether the algorithm will step back + // to a previous position as well as relax every time it + // detects a worsening of the metric; + // default is false: + bool m_BackTracking; + + // this is the number of successful iterations that must occur + // after which the algorithm will increase the step size; + // default is 1000000: + unsigned int m_PickUpPaceSteps; + + the_log_t * log_; + +private: + // disable the copy constructor and assignment operator: + RegularStepGradientDescentOptimizer2(const Self &); + void operator = (const Self &); +}; + +#if 0 +{ +#endif +} // end namespace itk + + +#endif // __itkRegularStepGradientDescentOptimizer2_h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index adde20b..15fb588 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -11,6 +11,11 @@ set(Nornir_SRCS IRTransaction.cxx IRLog.cxx IRThreadPool.cxx - ) + IRMosaicRefinementCommon.cxx + IRGridCommon.cxx + IRFFTCommon.cxx + IRGridTransform.cxx + itkRegularStepGradientDescentOptimizer2.cxx +) itk_module_add_library(Nornir ${Nornir_SRCS}) diff --git a/src/IRFFTCommon.cxx b/src/IRFFTCommon.cxx new file mode 100644 index 0000000..85cb62b --- /dev/null +++ b/src/IRFFTCommon.cxx @@ -0,0 +1,699 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : fft_common.cxx +// Author : Pavel A. Koshevoy +// Created : 2005/11/10 14:05 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : Helper functions for image alignment (registration) +// using phase correlation to find the translation vector. + +// local includes: +#include "IRFFTCommon.h" +#include "itkIRUtils.h" + +// system includes: +#include + +//---------------------------------------------------------------- +// DEBUG_COUNTER1 +// +#ifdef DEBUG_PDF +unsigned int DEBUG_COUNTER1 = 0; +unsigned int DEBUG_COUNTER2 = 0; +#endif + + +//---------------------------------------------------------------- +// find_maxima_cm +// +// The percentage below refers to the number of pixels that fall +// below the maxima. Thus, the number of pixels above the maxima +// is 1 - percentage. This way it is possible to specify a +// thresholding value without knowing anything about the image. +// +// The given image is thresholded, the resulting clusters/blobs +// are identified/classified, the center of mass of each cluster +// is treated as a maxima in the image. +// +unsigned int +find_maxima_cm(std::list & max_list, + const itk_image_t::Pointer & image, + const double percentage, + const the_text_t & prefix, + const the_text_t & suffix) +{ + typedef itk::ImageRegionConstIterator iter_t; + typedef itk::ImageRegionConstIteratorWithIndex itex_t; + typedef itk_image_t::IndexType index_t; + typedef std::list cluster_t; + + // local copy of the image that will be destroyed in the process: + itk_image_t::Pointer peaks = cast(image); + + // save(peaks, "Peaks.png"); + + // FIXME: +#ifdef DEBUG_PDF + if (prefix.size() != 0) + { + save(cast(remap_min_max(peaks, 0.0, 255.0)), + prefix + the_text_t("PDF") + suffix); + } +#endif + + // first find minimax/maxima of the image: + double v_min = std::numeric_limits::max(); + double v_max = -v_min; + + iter_t iter(peaks, peaks->GetLargestPossibleRegion()); + for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter) + { + const double v = iter.Get(); + v_min = std::min(v_min, v); + v_max = std::max(v_max, v); + } + + // calculate the min/max range: + const double v_rng = v_max - v_min; + + // NaN is the only number which is not equal to itself + if (v_rng == 0.0 || v_rng != v_rng || v_rng == std::numeric_limits::infinity()) + { + // there are no peaks in this image: + return 0; + } + + // build a histogram: + const unsigned int bins = 4096; + unsigned int pdf[bins] = { 0 }; + + for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter) + { + const double v = iter.Get(); + const unsigned int bin = (unsigned int)(double((v - v_min) / v_rng) * double(bins - 1)); + pdf[bin]++; + } + + // build the cumulative histogram: + unsigned int cdf[bins]; + cdf[0] = pdf[0]; + for (unsigned int i = 1; i < bins; i++) + { + cdf[i] = cdf[i - 1] + pdf[i]; + } + + // shortcuts: + itk_image_t::SizeType size = peaks->GetLargestPossibleRegion().GetSize(); + const unsigned int & w = size[0]; + const unsigned int & h = size[1]; + const double wh = double(w * h); + + // find the CDF bin that contains a given percentage of the total image: + double clip_min = 0.0; + for (unsigned int i = 1; i < bins; i++) + { + clip_min = v_min + (double(i) / double(bins - 1)) * v_rng; + if (double(cdf[i]) >= percentage * wh) + break; + } + + // threshold the peaks: + double background = clip_min - v_rng * 1e-3; + peaks = threshold(peaks, clip_min, v_max, background, v_max); + peaks = remap_min_max(peaks, 0.0, 1.0); + background = 0.0; + + // FIXME: +#ifdef DEBUG_CLUSTERS + if (prefix.size() != 0) + { + save(cast(remap_min_max(peaks, 0.0, 255.0)), + prefix + the_text_t("clusters") + suffix); + } +#endif + + // classify the clusters: + static const int stencil[][2] = { // 4 connected: + { 0, -1 }, + { -1, 0 }, + { 0, 1 }, + { 1, 0 }, + + // 8 connected: + { -1, -1 }, + { 1, 1 }, + { -1, 1 }, + { 1, -1 } + }; + + the_dynamic_array_t clusters; + the_dynamic_array_t bboxes; + std::vector cluster_map(w * h); + cluster_map.assign(w * h, ~0); + + itex_t itex(peaks, peaks->GetLargestPossibleRegion()); + for (itex.GoToBegin(); !itex.IsAtEnd(); ++itex) + { + const double v = itex.Get(); + + // skip over the background: + if (v <= background) + continue; + + const index_t index = itex.GetIndex(); + const unsigned int x = index[0]; + const unsigned int y = index[1]; + + // iterate over the neighborhood, collect the blob ids + // of the neighbors: + // James A: Changed to vector instead of a list... much faster + std::vector neighbors; + neighbors.reserve(8); + + for (unsigned int k = 0; k < 8; k++) + { + int u = x + stencil[k][0]; + int v = y + stencil[k][1]; + if ((unsigned int)(u) >= w || (unsigned int)(v) >= h) + continue; + + unsigned int cluster_id = cluster_map[u * h + v]; + if (cluster_id != (unsigned int)(~0)) + { + push_back_unique(neighbors, cluster_id); + } + } + + if (neighbors.empty()) + { + // make a new cluster: + clusters.append(cluster_t()); + bboxes.append(cluster_bbox_t()); + + unsigned int id = clusters.end_index(true); + cluster_map[x * h + y] = id; + clusters[id].push_back(index); + bboxes[id].update(x, y); + } + else + { + // add this pixel to the cluster: + unsigned int id = *(neighbors.begin()); + cluster_map[x * h + y] = id; + clusters[id].push_back(index); + bboxes[id].update(x, y); + + if (neighbors.size() > 1) + { + // merge the clusters into one (the first one): + // James A: Changed to vector instead of a list... much faster + std::vector::iterator bi = ++(neighbors.begin()); + for (; bi != neighbors.end(); ++bi) + { + unsigned int old_id = *bi; + bboxes[old_id].reset(); + + while (!clusters[old_id].empty()) + { + index_t ij = remove_head(clusters[old_id]); + + cluster_map[ij[0] * h + ij[1]] = id; + clusters[id].push_back(ij); + bboxes[id].update(ij[0], ij[1]); + } + } + } + } + } + + // merge the clusters that are broken up across the periodic boundary: + for (unsigned int i = 0; i < clusters.size(); i++) + { + cluster_t & cluster = clusters[i]; + if (cluster.empty()) + continue; + + for (std::list::iterator j = cluster.begin(); j != cluster.end(); ++j) + { + const index_t index = *j; + unsigned int x = (index[0] + w) % w; + unsigned int y = (index[1] + h) % h; + + for (unsigned int k = 0; k < 8; k++) + { + int dx = stencil[k][0]; + int dy = stencil[k][1]; + + // adjust for periodicity: + int u = (x + dx + w) % w; + int v = (y + dy + h) % h; + + unsigned int cluster_id = cluster_map[u * h + v]; + if (cluster_id == i || cluster_id == (unsigned int)(~0)) + continue; + + // figure out which boundaries this cluster was broken accross: + cluster_bbox_t & ba = bboxes[i]; + cluster_bbox_t & bb = bboxes[cluster_id]; + + bool merge_x = ((bb.max_[0] - ba.min_[0] > int(w / 2)) || (ba.max_[0] - bb.min_[0] > int(w / 2))); + + bool merge_y = ((bb.max_[1] - ba.min_[1] > int(h / 2)) || (ba.max_[1] - bb.min_[1] > int(h / 2))); + + int shift_x = (!merge_x) ? 0 : (ba.min_[0] <= 0) ? -w : w; + int shift_y = (!merge_y) ? 0 : (ba.min_[1] <= 0) ? -h : h; + + // merge the clusters into one (the first one): + cluster_t & neighbor = clusters[cluster_id]; + + bb.reset(); + while (!neighbor.empty()) + { + index_t ij = remove_head(neighbor); + cluster_map[ij[0] * h + ij[1]] = i; + + ij[0] += shift_x; + ij[1] += shift_y; + clusters[i].push_back(ij); + ba.update(ij[0], ij[1]); + } + } + } + } + + // FIXME: +#ifdef DEBUG_MARKERS + itk_image_t::Pointer markers = make_image(size, background); +#endif + + // calculate the center of mass for each cluster: + unsigned int num_peaks = 0; + for (unsigned int i = 0; i < clusters.size(); i++) + { + const cluster_t & cluster = clusters[i]; + if (cluster.empty()) + continue; + + double mx = 0.0; + double my = 0.0; + double mt = 0.0; + + for (std::list::const_iterator j = cluster.begin(); j != cluster.end(); ++j) + { + index_t ij = *j; + double x = double(ij[0]); + double y = double(ij[1]); + + // adjust index for periodicity: + if (x < 0) + ij[0] += w; + if (x >= w) + ij[0] -= w; + if (y < 0) + ij[1] += h; + if (y >= h) + ij[1] -= h; + + double m = peaks->GetPixel(ij); + mx += m * x; + my += m * y; + mt += m; + } + + double cm_x = mx / mt; + double cm_y = my / mt; + double m = mt / double(cluster.size()); + + // FIXME: +#ifdef DEBUG_MARKERS + mark(markers, pnt2d(cm_x, cm_y), m, 2, '+'); +#endif + + max_list.push_back(local_max_t(m, cm_x, cm_y, cluster.size())); + num_peaks++; + } + + // FIXME: +#ifdef DEBUG_MARKERS + save(cast(remap_min_max(markers, 0.0, 255.0)), + prefix + the_text_t("markings") + suffix, + false); +#endif + + // sort the max points so that the best candidate is first: + max_list.sort(std::greater()); + + return num_peaks; +} + + +//---------------------------------------------------------------- +// threshold_maxima +// +// Discard maxima whose mass is below a given threshold ratio +// of the total mass of all maxima: +// +void +threshold_maxima(std::list & max_list, const double threshold) +{ + double total_mass = 0.0; + for (std::list::iterator i = max_list.begin(); i != max_list.end(); ++i) + { + double mass = double((*i).area_) * ((*i).value_); + total_mass += mass; + } + + std::list new_list; + double threshold_mass = threshold * total_mass; + for (std::list::iterator i = max_list.begin(); i != max_list.end(); ++i) + { + double mass = double((*i).area_) * (*i).value_; + if (mass < threshold_mass) + continue; + + new_list.push_back(*i); + } + + max_list = new_list; +} + + +//---------------------------------------------------------------- +// reject_negligible_maxima +// +// Discard maxima that are worse than the best maxima by a factor +// greater than the given threshold ratio: +// +unsigned int +reject_negligible_maxima(std::list & max_list, const double threshold) +{ + double best_mass = 0.0; + for (std::list::iterator i = max_list.begin(); i != max_list.end(); ++i) + { + double mass = (*i).value_; + best_mass = std::max(best_mass, mass); + } + + unsigned int new_size = 0; + std::list new_list; + while (!max_list.empty()) + { + local_max_t lm = remove_head(max_list); + double mass = lm.value_; + if (best_mass / mass > threshold) + continue; + + new_list.push_back(lm); + new_size++; + } + + max_list.splice(max_list.end(), new_list); + return new_size; +} + + +//---------------------------------------------------------------- +// reject_negligible_overlap +// +void +reject_negligible_overlap(std::list & ol, const double threshold) +{ + double best_overlap = 0.0; + for (std::list::iterator i = ol.begin(); i != ol.end(); ++i) + { + double overlap = (*i).overlap_; + best_overlap = std::max(best_overlap, overlap); + } + + std::list new_list; + for (std::list::iterator i = ol.begin(); i != ol.end(); ++i) + { + double overlap = (*i).overlap_; + if (overlap == 0.0) + continue; + if (best_overlap / overlap > threshold) + continue; + + new_list.push_back(*i); + } + + ol = new_list; +} + + +//---------------------------------------------------------------- +// find_correlation +// +template <> +unsigned int +find_correlation(std::list & max_list, + const itk_image_t * fi, + const itk_image_t * mi, + + // low pass filter parameters + // (resampled data requires less smoothing): + double lp_filter_r, + double lp_filter_s, + const double overlap_min, + const double overlap_max) +{ + itk_image_t::SizeType max_sz = calc_padding(fi, mi); + + typedef itk_image_t::SizeType sz_t; + typedef itk_image_t::RegionType rn_t; + + rn_t fi_region = fi->GetLargestPossibleRegion(); + sz_t fi_size = fi_region.GetSize(); + + rn_t mi_region = mi->GetLargestPossibleRegion(); + sz_t mi_size = mi_region.GetSize(); + + itk_image_t::ConstPointer z0 = + (fi_size[0] == max_sz[0] && fi_size[1] == max_sz[1]) ? fi : pad(fi, max_sz); + itk_image_t::ConstPointer z1 = + (mi_size[0] == max_sz[0] && mi_size[1] == max_sz[1]) ? mi : pad(mi, max_sz); + + fft_data_t f0; + fft(z0, f0); + f0.apply_lp_filter(lp_filter_r, lp_filter_s); + + fft_data_t f1; + fft(z1, f1); + f1.apply_lp_filter(lp_filter_r, lp_filter_s); + + const unsigned int & nx = f0.nx(); + const unsigned int & ny = f0.ny(); + fft_data_t P(nx, ny); + + // Blank out areas outside the min or max overlap + unsigned int xLeftBorderStart = std::ceil(nx * overlap_min); + unsigned int xRightBorderStart = nx - xLeftBorderStart; + unsigned int yLowBorderStart = std::ceil(ny * overlap_min); + unsigned int yHighBorderStart = ny - yLowBorderStart; + + + unsigned int xLeftBorderCutoff = std::ceil(nx * overlap_max); + unsigned int xRightBorderCutoff = nx - xLeftBorderCutoff; + unsigned int yLowBorderCutoff = std::ceil(ny * overlap_max); + unsigned int yHighBorderCutoff = ny - yLowBorderCutoff; + + for (unsigned int x = 0; x < nx; x++) + { + for (unsigned int y = 0; y < ny; y++) + { + // TODO: Skip pixels we know can't possibly be overlapping + /* + if(!(y <= yLowBorderCutoff || + y >= yHighBorderCutoff || + x <= xLeftBorderCutoff || + x >= xRightBorderCutoff)) + { + P(x, y) = 0; + continue; + } + */ +#if 1 + // Girod-Kuo, normalized cross power spectrum, + // corresponds to phase correlation in spatial domain: + fft_complex_t p10 = f1(x, y) * std::conj(f0(x, y)); + P(x, y) = _div(p10, _add(std::sqrt(p10 * std::conj(p10)), 1e-8f)); +#else + // cross power spectrum, + // corresponds to cross correlation in spatial domain: + P(x, y) = f1(x, y) * std::conj(f0(x, y)); +#endif + } + } + + // resampled data produces less noisy PDF and requires less smoothing: + P.apply_lp_filter(lp_filter_r * 0.8, lp_filter_s); + + // calculate the displacement probability density function: + fft_data_t ifft_P; + +#ifndef NDEBUG // get around an annoying compiler warning: + bool ok = +#endif + ifft(P, ifft_P); + assert(ok); + + itk_image_t::Pointer PDF = ifft_P.real(); + + itk_image_t::PixelType min; + itk_image_t::PixelType max; + image_min_max(PDF.GetPointer(), min, max); + + // save(cast(remap_min_max(PDF, 0.0, 255.0)), + // "Prefill.tif"); + + // Set areas that cannot be a match to zero + // TODO: //Need to fill with an ellipse to get mins/maxs correct + + // if(xLeftBorderCutoff < xRightBorderCutoff && + // yLowBorderCutoff < yHighBorderCutoff) + // { + // fill(PDF, xLeftBorderCutoff, yLowBorderCutoff, xRightBorderCutoff-xLeftBorderCutoff, + //yHighBorderCutoff-yLowBorderCutoff, min); + // } + // + // #ifdef DEBUG + // save(cast(remap_min_max(PDF, 0.0, 255.0)), + //"Prefill.tif"); + // #endif + + int PixelsInOverlapZone = 0; + itk_image_t::IndexType iPixel; + for (iPixel[1] = 0; iPixel[1] <= ny / 2; iPixel[1]++) + { + vec2d_t pt; + for (iPixel[0] = 0; iPixel[0] <= nx / 2; iPixel[0]++) + { + pt[0] = iPixel[0]; + pt[1] = iPixel[1]; + + double overlap = OverlapPercent(fi_size, pt); + if (overlap >= overlap_min && overlap <= overlap_max) + { + PixelsInOverlapZone += 4; + continue; + } + + pt[0] = nx - iPixel[0]; + pt[1] = iPixel[1]; + + overlap = OverlapPercent(fi_size, pt); + if (overlap >= overlap_min && overlap <= overlap_max) + { + PixelsInOverlapZone += 4; + continue; + } + + pt[0] = iPixel[0]; + pt[1] = ny - iPixel[1]; + + overlap = OverlapPercent(fi_size, pt); + if (overlap >= overlap_min && overlap <= overlap_max) + { + PixelsInOverlapZone += 4; + continue; + } + + pt[0] = nx - iPixel[0]; + pt[1] = ny - iPixel[1]; + + overlap = OverlapPercent(fi_size, pt); + if (overlap >= overlap_min && overlap <= overlap_max) + { + PixelsInOverlapZone += 4; + continue; + } + + // If we got here the pixel can't be overlapping + itk_image_t::IndexType iSetPixel = iPixel; + PDF->SetPixel(iSetPixel, min); + + iSetPixel[0] = (nx - 1) - iPixel[0]; + PDF->SetPixel(iSetPixel, min); + + iSetPixel[0] = iPixel[0]; + iSetPixel[1] = (ny - 1) - iPixel[1]; + PDF->SetPixel(iSetPixel, min); + + iSetPixel[0] = (nx - 1) - iPixel[0]; + iSetPixel[1] = (ny - 1) - iPixel[1]; + PDF->SetPixel(iSetPixel, min); + /* + if(overlap < overlap_min) + PDF->SetPixel(iPixel, min); + + else if(overlap > overlap_max) + PDF->SetPixel(iPixel, min); + */ + } + } + + /* + itk_image_t::IndexType ix; + ix[0] = 0; + ix[1] = 0; + PDF->SetPixel(ix, min); + ix[0] = nx-1; + ix[1] = 0; + PDF->SetPixel(ix, min); + ix[0] = 0; + ix[1] = ny-1; + PDF->SetPixel(ix, min); + ix[0] = nx-1; + ix[1] = ny-1; + PDF->SetPixel(ix, min); + */ + // fill(PDF, 0, 0, xLeftBorderStart, ny, min); + // ill(PDF, xRightBorderStart, 0, nx-xRightBorderStart, ny, min); + // fill(PDF, xLeftBorderStart, 0, xRightBorderStart-xLeftBorderStart, yLowBorderStart, min); + // fill(PDF, xLeftBorderStart, yHighBorderStart, xRightBorderStart-xLeftBorderStart, + //ny-yHighBorderStart, min); + + // #ifdef DEBUG + // save(cast(remap_min_max(PDF, 0.0, 255.0)), + //"Postfill.tif"); + // #endif + + // look for the maxima in the PDF, + // TODO: Should we count the dead regions towards our histogram because they'll be examined? + // double area = double(max_sz[0] * max_sz[1]); + // if(xLeftBorderCutoff < xRightBorderCutoff && + // yLowBorderCutoff < yHighBorderCutoff) + // area -= ((xRightBorderCutoff-xLeftBorderCutoff) * (yHighBorderCutoff-yLowBorderCutoff)); + double area = PixelsInOverlapZone; + + // a minimum of 5 pixels and a maximum of 64 pixels may be attributed + // to local maxima in the image: + // double fraction = std::min(64.0 / area, std::max(5.0 / area, 1e-3)); + double fraction = std::min(64.0 / area, std::max(5.0 / area, 1e-2)); + + // the entire image should never be treated as a maxima cluster: + assert(fraction < 1.0); + + // find the maxima clusters: + return find_maxima_cm(max_list, PDF, 1.0 - fraction); +} diff --git a/src/IRGridCommon.cxx b/src/IRGridCommon.cxx new file mode 100644 index 0000000..8405bc1 --- /dev/null +++ b/src/IRGridCommon.cxx @@ -0,0 +1,333 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : grid_common.cxx +// Author : Pavel A. Koshevoy +// Created : Wed Jan 10 09:31:00 MDT 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : code used to refine mesh transform control points + +// the includes: +#include "IRGridCommon.h" + + +//---------------------------------------------------------------- +// setup_grid_transform +// +bool +setup_grid_transform(the_grid_transform_t & transform, + unsigned int rows, + unsigned int cols, + const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const mask_t * tile_mask, + base_transform_t::ConstPointer mosaic_to_tile, + unsigned int max_iterations, + double min_step_scale, + double min_error_sqrd, + unsigned int pick_up_pace_steps) +{ + const itk::GridTransform * gt = dynamic_cast(mosaic_to_tile.GetPointer()); + + std::vector xy_arr((rows + 1) * (cols + 1)); + std::vector xy_apx((rows + 1) * (cols + 1)); + + image_t::Pointer dx = make_image(cols + 1, rows + 1, 1.0, 0.0); + image_t::Pointer dy = make_image(cols + 1, rows + 1, 1.0, 0.0); + + typedef itk::LegendrePolynomialTransform approx_transform_t; + + // the mosaic to tile transform is typically more stable: + approx_transform_t::Pointer mosaic_to_tile_approx; + + if (gt == NULL) + { + mosaic_to_tile_approx = approx_transform(tile_min, + tile_max, + tile_mask, + mosaic_to_tile.GetPointer(), + 16, // samples per edge + 1, // landmark generator version + true); // iterative refinement + } + + // temporaries: + vec2d_t tile_ext = tile_max - tile_min; + + bool ReturnVal = true; + +#pragma omp parallel for + for (int row = 0; row <= (int)rows; row++) + { + pnt2d_t uv; + pnt2d_t pq; + + pq[1] = double(row) / double(rows); + uv[1] = tile_min[1] + tile_ext[1] * pq[1]; + for (int col = 0; col <= (int)cols; col++) + { + pq[0] = double(col) / double(cols); + uv[0] = tile_min[0] + tile_ext[0] * pq[0]; + + // shortcut: + unsigned int index = row * (cols + 1) + col; + pnt2d_t & xy = xy_arr[index]; + pnt2d_t & xy_approx = xy_apx[index]; + + if (gt == NULL) + { + // general transform: + if (!find_inverse(tile_min, + tile_max, + mosaic_to_tile_approx.GetPointer(), + uv, + xy_approx, + max_iterations, + min_step_scale, + min_error_sqrd, + pick_up_pace_steps)) + { + // we are screwed: + ReturnVal |= false; + break; + } + + if (!find_inverse(tile_min, + tile_max, + mosaic_to_tile.GetPointer(), + uv, + xy, + max_iterations, + min_step_scale, + min_error_sqrd, + pick_up_pace_steps)) + { + xy = xy_approx; + } + + pnt2d_t uv2 = mosaic_to_tile->TransformPoint(xy); + + // verify that the point maps back correctly within some tolerance: + vec2d_t e_uv = uv2 - uv; + + // verify that the approximate and exact aren't too far apart: + vec2d_t e_xy = xy_approx - xy; + + double e_uv_absolute = e_uv.GetSquaredNorm(); + + // FIXME: this is an idea -- if the exact transform give wildely + // differing result from the approximate transform, we may be able + // remove such outliers via a median filter. To do that, we have to + // maintain an image of the approximate/exact result differences: + image_t::IndexType ix; + ix[0] = col; + ix[1] = row; + dx->SetPixel(ix, e_xy[0]); + dy->SetPixel(ix, e_xy[1]); + + // FIXME: this is a temporary crutch, the method outlined above + // should be used instead: + static const double uv_tolerance = 1e-6; + if (e_uv_absolute > uv_tolerance) + { + xy = xy_approx; + } + } + else + { + // discontinuous transform -- this is a more stable way to resample + // the transformation mesh: + bool ok = gt->transform_.transform_inv(pq, xy); +#if !defined(__APPLE__) + assert(ok); + assert(xy[0] == xy[0] && xy[1] == xy[1]); + ReturnVal &= ok; +#endif + } + } + } + + if (!ReturnVal) + return false; + +#if 0 + save(cast + (remap_min_max(dx)), + "init-error-x.tif"); + + save(cast + (remap_min_max(dy)), + "init-error-y.tif"); +#endif + + transform.setup(rows, cols, tile_min, tile_max, xy_arr); + return true; +} + +//---------------------------------------------------------------- +// setup_mesh_transform +// +bool +setup_mesh_transform(the_mesh_transform_t & transform, + unsigned int rows, + unsigned int cols, + const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const mask_t * tile_mask, + base_transform_t::ConstPointer mosaic_to_tile, + unsigned int max_iterations, + double min_step_scale, + double min_error_sqrd, + unsigned int pick_up_pace_steps) +{ + const itk::GridTransform * gt = dynamic_cast(mosaic_to_tile.GetPointer()); + + std::vector xy_arr((rows + 1) * (cols + 1)); + std::vector xy_apx((rows + 1) * (cols + 1)); + + image_t::Pointer dx = make_image(cols + 1, rows + 1, 1.0, 0.0); + image_t::Pointer dy = make_image(cols + 1, rows + 1, 1.0, 0.0); + + typedef itk::LegendrePolynomialTransform approx_transform_t; + + // the mosaic to tile transform is typically more stable: + approx_transform_t::Pointer mosaic_to_tile_approx; + + if (gt == NULL) + { + mosaic_to_tile_approx = approx_transform(tile_min, + tile_max, + tile_mask, + mosaic_to_tile.GetPointer(), + 16, // samples per edge + 1, // landmark generator version + true); // iterative refinement + } + + // temporaries: + vec2d_t tile_ext = tile_max - tile_min; + pnt2d_t uv; + pnt2d_t pq; + + std::vector uv_list((rows + 1) * (cols + 1)); + + for (unsigned int row = 0; row <= rows; row++) + { + pq[1] = double(row) / double(rows); + uv[1] = tile_min[1] + tile_ext[1] * pq[1]; + for (unsigned int col = 0; col <= cols; col++) + { + pq[0] = double(col) / double(cols); + uv[0] = tile_min[0] + tile_ext[0] * pq[0]; + + // shortcut: + unsigned int index = row * (cols + 1) + col; + pnt2d_t & xy = xy_arr[index]; + pnt2d_t & xy_approx = xy_apx[index]; + uv_list[index] = pq; + + if (gt == NULL) + { + // general transform: + if (!find_inverse(tile_min, + tile_max, + mosaic_to_tile_approx.GetPointer(), + uv, + xy_approx, + max_iterations, + min_step_scale, + min_error_sqrd, + pick_up_pace_steps)) + { + // we are screwed: + return false; + } + + if (!find_inverse(tile_min, + tile_max, + mosaic_to_tile.GetPointer(), + uv, + xy, + max_iterations, + min_step_scale, + min_error_sqrd, + pick_up_pace_steps)) + { + xy = xy_approx; + } + + pnt2d_t uv2 = mosaic_to_tile->TransformPoint(xy); + + // verify that the point maps back correctly within some tolerance: + vec2d_t e_uv = uv2 - uv; + + // verify that the approximate and exact aren't too far apart: + vec2d_t e_xy = xy_approx - xy; + + double e_uv_absolute = e_uv.GetSquaredNorm(); + + // FIXME: this is an idea -- if the exact transform give wildely + // differing result from the approximate transform, we may be able + // remove such outliers via a median filter. To do that, we have to + // maintain an image of the approximate/exact result differences: + image_t::IndexType ix; + ix[0] = col; + ix[1] = row; + dx->SetPixel(ix, e_xy[0]); + dy->SetPixel(ix, e_xy[1]); + + // FIXME: this is a temporary crutch, the method outlined above + // should be used instead: + static const double uv_tolerance = 1e-6; + if (e_uv_absolute > uv_tolerance) + { + xy = xy_approx; + } + } + else + { + // discontinuous transform -- this is a more stable way to resample + // the transformation mesh: + bool ok = gt->transform_.transform_inv(pq, xy); +#if !defined(__APPLE__) + assert(ok); + assert(xy[0] == xy[0] && xy[1] == xy[1]); +#endif + if (!ok) + return false; + } + } + } + +#if 0 + save(cast + (remap_min_max(dx)), + "init-error-x.tif"); + + save(cast + (remap_min_max(dy)), + "init-error-y.tif"); +#endif + + transform.setup(tile_min, tile_max, uv_list, xy_arr); + return true; +} diff --git a/src/IRGridTransform.cxx b/src/IRGridTransform.cxx new file mode 100644 index 0000000..96f3279 --- /dev/null +++ b/src/IRGridTransform.cxx @@ -0,0 +1,1487 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: nil -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_grid_transform.cxx +// Author : Pavel A. Koshevoy +// Created : Thu Nov 30 13:37:18 MST 2006 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A discontinuous transform -- a uniform grid of vertices is +// mapped to an image. At each vertex, in addition to image +// space coordinates, a second set of coordinates is stored. +// This is similar to texture mapped OpenGL triangle meshes, +// where the texture coordinates correspond to the image space +// vertex coordinates. + +// local includes: +#include "IRGridTransform.h" + +// system includes: +#include +#include +#include + + +//---------------------------------------------------------------- +// GRID_DENSITY +// +static const int GRID_DENSITY = 1; + +//---------------------------------------------------------------- +// EPSILON +// +static const double EPSILON = 1e-6; + + +//---------------------------------------------------------------- +// WARPING +// +#ifdef DEBUG_WARPING +bool WARPING = false; +#endif + +//---------------------------------------------------------------- +// triangle_t::triangle_t +// +triangle_t::triangle_t() +{ + vertex_[0] = ~0; + vertex_[1] = ~0; + vertex_[2] = ~0; +} + +//---------------------------------------------------------------- +// triangle_intersect +// +inline static bool +triangle_intersect(const double * pwb, + const double * pwc, + const double & pt_x, + const double & pt_y, + double * barycentric_coords) +{ +#define wa barycentric_coords[0] +#define wb barycentric_coords[1] +#define wc barycentric_coords[2] + + wb = pwb[0] * pt_x + pwb[1] * pt_y + pwb[2]; + wc = pwc[0] * pt_x + pwc[1] * pt_y + pwc[2]; + + // corner cases -- clamp when within tolerance: + if (wb < 0 && wb > -EPSILON) + wb = 0; + if (wb > 1 && (wb - 1) < EPSILON) + wb = 1; + if (wc < 0 && wc > -EPSILON) + wc = 0; + if (wc > 1 && (wc - 1) < EPSILON) + wc = 1; + + wa = 1.0 - (wb + wc); + if (wa < 0 && wa > -EPSILON) + wa = 0; + if (wa > 1 && (wa - 1) < EPSILON) + wa = 1; + + return (wa >= 0.0 && wb >= 0.0 && wc >= 0.0); + +#undef wa +#undef wb +#undef wc +} + +//---------------------------------------------------------------- +// triangle_t::intersect +// +bool +triangle_t::xy_intersect(const vertex_t * v_arr, const pnt2d_t & xy, pnt2d_t & uv) const +{ +#define wa barycentric_coords[0] +#define wb barycentric_coords[1] +#define wc barycentric_coords[2] + + double barycentric_coords[3]; + if (!triangle_intersect(xy_pwb, xy_pwc, xy[0], xy[1], barycentric_coords)) + { + return false; + } + + const pnt2d_t & A = v_arr[vertex_[0]].uv_; + const pnt2d_t & B = v_arr[vertex_[1]].uv_; + const pnt2d_t & C = v_arr[vertex_[2]].uv_; + + uv[0] = A[0] * wa + B[0] * wb + C[0] * wc; + uv[1] = A[1] * wa + B[1] * wb + C[1] * wc; + +#ifdef DEBUG_WARPING + static const double THRESHOLD = 5e-2; + if (WARPING && (fabs(wa) < THRESHOLD || fabs(wb) < THRESHOLD || fabs(wc) < THRESHOLD)) + { + uv[0] = 0.0; + uv[1] = 0.0; + } +#endif // DEBUG_WARPING + + return true; + +#undef wa +#undef wb +#undef wc +} + +//---------------------------------------------------------------- +// triangle_t::uv_intersect +// +bool +triangle_t::uv_intersect(const vertex_t * v_arr, const pnt2d_t & uv, pnt2d_t & xy) const +{ +#define wa barycentric_coords[0] +#define wb barycentric_coords[1] +#define wc barycentric_coords[2] + + double barycentric_coords[3]; + if (!triangle_intersect(uv_pwb, uv_pwc, uv[0], uv[1], barycentric_coords)) + { +#if 0 + // for debugging: + static int count = 0; + count++; + printf("%3i. %+e %+e %+e, %e\n", count, wa, wb, wc, + fabs(wa) + fabs(wb) + fabs(wc)); + triangle_intersect(uv_pwb, uv_pwc, uv[0], uv[1], barycentric_coords); +#endif + + return false; + } + + const pnt2d_t & A = v_arr[vertex_[0]].xy_; + const pnt2d_t & B = v_arr[vertex_[1]].xy_; + const pnt2d_t & C = v_arr[vertex_[2]].xy_; + + xy[0] = A[0] * wa + B[0] * wb + C[0] * wc; + xy[1] = A[1] * wa + B[1] * wb + C[1] * wc; + + return true; + +#undef wa +#undef wb +#undef wc +} + + +//---------------------------------------------------------------- +// the_acceleration_grid_t::the_acceleration_grid_t +// +the_acceleration_grid_t::the_acceleration_grid_t() + : rows_(0) + , cols_(0) +{ + // reset the grid bounding box: + xy_min_[0] = std::numeric_limits::max(); + xy_min_[1] = xy_min_[0]; + xy_ext_[0] = 0.0; + xy_ext_[1] = 0.0; +} + +//---------------------------------------------------------------- +// the_acceleration_grid_t::xy_cell +// +unsigned int +the_acceleration_grid_t::xy_cell(const pnt2d_t & xy) const +{ + // find where in the grid the point lands: + double a = (xy[0] - xy_min_[0]) / xy_ext_[0]; + double b = (xy[1] - xy_min_[1]) / xy_ext_[1]; + if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0) + { + double c = a * double(cols_); + double r = b * double(rows_); + + unsigned int col = c >= cols_ ? cols_ - 1 : c; + unsigned int row = r >= rows_ ? rows_ - 1 : r; + return row * cols_ + col; + } + + return ~0; +} + +//---------------------------------------------------------------- +// the_acceleration_grid_t::xy_triangle +// +unsigned int +the_acceleration_grid_t::xy_triangle(const pnt2d_t & xy, pnt2d_t & uv) const +{ + unsigned int cell_id = xy_cell(xy); + + if (cell_id != (unsigned int)(~0)) + { + // shortcuts: + const std::vector & cell = xy_[cell_id]; + const vertex_t * v_arr = &(mesh_[0]); + + // check each candidate triangle for intersection: + for (int iter = 0; iter < (int)cell.size(); iter++) + { + unsigned int iTri = cell[iter]; + const triangle_t & tri = tri_[iTri]; + if (tri.xy_intersect(v_arr, xy, uv)) + { + return iTri; + } + } + } + + return ~0; +} + +//---------------------------------------------------------------- +// the_acceleration_grid_t::uv_cell +// +unsigned int +the_acceleration_grid_t::uv_cell(const pnt2d_t & uv) const +{ + // find where in the grid the point lands: + double a = uv[0]; + double b = uv[1]; + if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0) + { + double c = std::min(double(cols_ - 1), a * double(cols_)); + double r = std::min(double(rows_ - 1), b * double(rows_)); + + unsigned int col = (unsigned int)(floor(c)); + unsigned int row = (unsigned int)(floor(r)); + + return row * cols_ + col; + } + + return ~0; +} + +//---------------------------------------------------------------- +// the_acceleration_grid_t::uv_triangle +// +unsigned int +the_acceleration_grid_t::uv_triangle(const pnt2d_t & uv, pnt2d_t & xy) const +{ + unsigned int cell_id = uv_cell(uv); + + if (cell_id != (unsigned int)(~0)) + { + // shortcuts: + const std::vector & cell = uv_[cell_id]; + const vertex_t * v_arr = &(mesh_[0]); + + // check each candidate triangle for intersection: + for (unsigned int iter = 0; iter < cell.size(); iter++) + { + unsigned int iTri = cell[iter]; + const triangle_t & tri = tri_[iTri]; + if (tri.uv_intersect(v_arr, uv, xy)) + { + return iTri; + } + } + } + + return ~0; +} + +//---------------------------------------------------------------- +// the_acceleration_grid_t::update +// +void +the_acceleration_grid_t::update(const vec2d_t * xy_shift) +{ + unsigned int num_verts = mesh_.size(); + vertex_t * v_arr = &(mesh_[0]); + for (unsigned int i = 0; i < num_verts; i++) + { + // cerr << v_arr[i].xy_ << " -> "; + v_arr[i].xy_ += xy_shift[i]; + // cerr << v_arr[i].xy_ << endl; + } + + rebuild(); +} + +//---------------------------------------------------------------- +// the_acceleration_grid_t::shift +// +void +the_acceleration_grid_t::shift(const vec2d_t & xy_shift) +{ + unsigned int num_verts = mesh_.size(); + vertex_t * v_arr = &(mesh_[0]); + for (unsigned int i = 0; i < num_verts; i++) + { + // cerr << v_arr[i].xy_ << " -> "; + v_arr[i].xy_ += xy_shift; + // cerr << v_arr[i].xy_ << endl; + } + + rebuild(); +} + +//---------------------------------------------------------------- +// the_acceleration_grid_t::resize +// +void +the_acceleration_grid_t::resize(unsigned int rows, unsigned int cols) +{ + // reset the grid bounding box: + xy_min_[0] = std::numeric_limits::max(); + xy_min_[1] = xy_min_[0]; + xy_ext_[0] = 0.0; + xy_ext_[1] = 0.0; + + rows_ = rows; + cols_ = cols; + xy_.resize(rows_ * cols_); + uv_.resize(rows_ * cols_); +} + +//---------------------------------------------------------------- +// the_acceleration_grid_t::rebuild +// +void +the_acceleration_grid_t::rebuild() +{ + // reset the grid bounding box and destroy the the old grid: + xy_min_[0] = std::numeric_limits::max(); + xy_min_[1] = xy_min_[0]; + + pnt2d_t xy_max; + xy_max[0] = -xy_min_[0]; + xy_max[1] = -xy_min_[1]; + + // find the new grid bounding box: + unsigned int num_verts = mesh_.size(); + vertex_t * v_arr = &(mesh_[0]); + for (unsigned int i = 0; i < num_verts; i++) + { + update_bbox(xy_min_, xy_max, v_arr[i].xy_); + } + xy_ext_ = xy_max - xy_min_; + +#if 0 + cerr << "xy_min: " << xy_min_ << endl + << "xy_max: " << xy_max << endl + << "xy_ext: " << xy_ext_ << endl; +#endif + + // reset to empty grid: + xy_.assign(xy_.size(), std::vector()); + uv_.assign(uv_.size(), std::vector()); + + xy_.reserve(10); + uv_.reserve(10); + + // add triangles to the grid: + unsigned int num_triangles = tri_.size(); + for (unsigned int i = 0; i < num_triangles; i++) + { + update_grid(i); + } + +#if 0 + int count = 0; + for (unsigned int row = 0; row < rows_; row++) + { + for (unsigned int col = 0; col < cols_; col++) + { + count += xy_[row * cols_ + col].size(); + cout << setw(3) << xy_[row * cols_ + col].size(); + } + cout << endl; + } + + assert(count > 0); +#endif +} + +//---------------------------------------------------------------- +// intersect_lines_2d +// +// Find the intersection parameters (u, v) of 2 parametric lines +// A = a + b * u +// B = c + d * v +// +static bool +intersect_lines_2d( // start point A: + double ax, + double ay, + + // direction vector A: + double bx, + double by, + + // start point B: + double cx, + double cy, + + // direction vector B: + double dx, + double dy, + + // results: + double & u, + double & v) +{ + // solve the linear system Ax = b: + // + // [bx -dx] [u] = [cx - ax] + // [by -dy] [v] [cy - ay] + // + + double det_A = dx * by - bx * dy; + if (det_A == 0.0) + return false; + + double det_A_inv = 1.0 / det_A; + double A_inv[][2] = { { -dy * det_A_inv, dx * det_A_inv }, { -by * det_A_inv, bx * det_A_inv } }; + + double b[] = { cx - ax, cy - ay }; + u = A_inv[0][0] * b[0] + A_inv[0][1] * b[1]; + v = A_inv[1][0] * b[0] + A_inv[1][1] * b[1]; + + return true; +} + +//---------------------------------------------------------------- +// intersect_bbox_triangle +// +static bool +intersect_bbox_triangle(const pnt2d_t & min, + const pnt2d_t & max, + const pnt2d_t & A, + const pnt2d_t & B, + const pnt2d_t & C) +{ + double lines_a[][4] = { // horizontal: + { min[0], min[1], max[0] - min[0], 0.0 }, + { min[0], max[1], max[0] - min[0], 0.0 }, + // vertical: + { min[0], min[1], 0.0, max[1] - min[1] }, + { max[0], min[1], 0.0, max[1] - min[1] } + }; + + double lines_b[][4] = { { A[0], A[1], B[0] - A[0], B[1] - A[1] }, + { A[0], A[1], C[0] - A[0], C[1] - A[1] }, + { B[0], B[1], C[0] - B[0], C[1] - B[1] } }; + + double u; + double v; + for (unsigned int i = 0; i < 4; i++) + { + for (unsigned int j = 0; j < 3; j++) + { + if (intersect_lines_2d(lines_a[i][0], + lines_a[i][1], + lines_a[i][2], + lines_a[i][3], + + lines_b[j][0], + lines_b[j][1], + lines_b[j][2], + lines_b[j][3], + + u, + v)) + { + if (u >= 0.0 && u <= 1.0 && v >= 0.0 && v <= 1.0) + { + return true; + } + } + } + } + + return false; +} + +//---------------------------------------------------------------- +// update_grid +// +static void +update_grid( // the acceleration grid: + std::vector * grid, + const unsigned int rows, + const unsigned int cols, + + // acceleration grid bounding box: + const pnt2d_t & grid_min, + const vec2d_t & grid_ext, + + // the triangle being added to the grid: + const pnt2d_t & A, + const pnt2d_t & B, + const pnt2d_t & C, + const unsigned int tri_id, + + // calculate fast barycentric coordinate calculation coefficients: + double * pwb, + double * pwc) +{ + // precompute the barycentric calculation coefficients: + { + vec2d_t b = B - A; + vec2d_t c = C - A; + double bycx_bxcy = b[1] * c[0] - b[0] * c[1]; + + pwb[0] = -c[1] / bycx_bxcy; + pwb[1] = c[0] / bycx_bxcy; + pwb[2] = (c[1] * A[0] - c[0] * A[1]) / bycx_bxcy; + + pwc[0] = b[1] / bycx_bxcy; + pwc[1] = -b[0] / bycx_bxcy; + pwc[2] = (b[0] * A[1] - b[1] * A[0]) / bycx_bxcy; + } + + pnt2d_t min = A; + pnt2d_t max = min; + update_bbox(min, max, B); + update_bbox(min, max, C); + + const double & gw = grid_ext[0]; + double a[] = { (min[0] - grid_min[0]) / gw, (max[0] - grid_min[0]) / gw }; + + unsigned int c[] = { std::min(cols - 1, (unsigned int)(floor(a[0] * double(cols)))), + std::min(cols - 1, (unsigned int)(floor(a[1] * double(cols)))) }; + + const double & gh = grid_ext[1]; + double b[] = { (min[1] - grid_min[1]) / gh, (max[1] - grid_min[1]) / gh }; + + unsigned int r[] = { std::min(rows - 1, (unsigned int)(floor(b[0] * double(rows)))), + std::min(rows - 1, (unsigned int)(floor(b[1] * double(rows)))) }; + + // temporary barycentric coordinates of point inside the triangle: + double barycentric_coords[3]; + + for (unsigned int row = r[0]; row <= r[1]; row++) + { + for (unsigned int col = c[0]; col <= c[1]; col++) + { + unsigned int i = row * cols + col; + std::vector & cell = grid[i]; + + // calculate the bounding box of this grid cell: + pnt2d_t min(grid_min); + min[0] += gw * double(col) / double(cols); + min[1] += gh * double(row) / double(rows); + + pnt2d_t max(grid_min); + max[0] += gw * double(col + 1) / double(cols); + max[1] += gh * double(row + 1) / double(rows); + + if (inside_bbox(min, max, A) || inside_bbox(min, max, B) || inside_bbox(min, max, C) || + triangle_intersect(pwb, pwc, min[0], min[1], barycentric_coords) || + triangle_intersect(pwb, pwc, max[0], min[1], barycentric_coords) || + triangle_intersect(pwb, pwc, max[0], max[1], barycentric_coords) || + triangle_intersect(pwb, pwc, min[0], max[1], barycentric_coords) || + intersect_bbox_triangle(min, max, A, B, C)) + { + // James A: Modifying this to be a bit more selective about which triangles are added to which grid cells. If a + // triangle only shares an edge with the top or right cell wall we do not add it + unsigned int numVertsInsideBBox = 0; + numVertsInsideBBox += inside_bbox(min, max, A); + numVertsInsideBBox += inside_bbox(min, max, B); + numVertsInsideBBox += inside_bbox(min, max, C); + + // If all three verts are inside the box then this it is more likely this triangle will match a test than a + // triangle with only one or two verts inside the box. Put it at the front of the line + if (numVertsInsideBBox == 3) + cell.insert(cell.begin(), tri_id); + else + cell.push_back(tri_id); + } + } + } +} + +//---------------------------------------------------------------- +// the_acceleration_grid_t::update_grid +// +void +the_acceleration_grid_t::update_grid(unsigned int t_idx) +{ + // shortcuts: + triangle_t & tri = tri_[t_idx]; + const vertex_t * v_arr = &(mesh_[0]); + const vertex_t & v0 = v_arr[tri.vertex_[0]]; + const vertex_t & v1 = v_arr[tri.vertex_[1]]; + const vertex_t & v2 = v_arr[tri.vertex_[2]]; + + // update the xy-grid: + ::update_grid(&(xy_[0]), + rows_, + cols_, + + // xy-grid bbox: + xy_min_, + xy_ext_, + + // xy-triangle: + v0.xy_, + v1.xy_, + v2.xy_, + + t_idx, + tri.xy_pwb, + tri.xy_pwc); + + // update the uv-grid: + ::update_grid(&(uv_[0]), + rows_, + cols_, + + // uv-grid bbox: + pnt2d(0, 0), + vec2d(1, 1), + + // xy-triangle: + v0.uv_, + v1.uv_, + v2.uv_, + + t_idx, + tri.uv_pwb, + tri.uv_pwc); +} + + +//---------------------------------------------------------------- +// the_base_triangle_transform_t::transform +// +bool +the_base_triangle_transform_t::transform(const pnt2d_t & xy, pnt2d_t & uv) const +{ + unsigned int t_id = grid_.xy_triangle(xy, uv); + + if (t_id == (unsigned int)(~0)) + { + uv[0] = std::numeric_limits::quiet_NaN(); + uv[1] = uv[0]; + } + + return t_id != (unsigned int)(~0); +} + +//---------------------------------------------------------------- +// the_base_triangle_transform_t::transform_inv +// +bool +the_base_triangle_transform_t::transform_inv(const pnt2d_t & uv, pnt2d_t & xy) const +{ + unsigned int t_id = grid_.uv_triangle(uv, xy); + + if (t_id == (unsigned int)(~0)) + { +#if 0 + // for debugging: + grid_.uv_triangle(uv, xy); +#endif + + xy[0] = std::numeric_limits::quiet_NaN(); + xy[1] = xy[0]; + } + + return t_id != (unsigned int)(~0); +} + +//---------------------------------------------------------------- +// the_base_triangle_transform_t::jacobian +// +bool +the_base_triangle_transform_t::jacobian(const pnt2d_t & P, unsigned int * idx, double * jac) const +{ + pnt2d_t uv; + unsigned int t_id = grid_.xy_triangle(P, uv); + if (t_id == (unsigned int)(~0)) + return false; + + const triangle_t & tri = grid_.tri_[t_id]; + idx[0] = tri.vertex_[0]; + idx[1] = tri.vertex_[1]; + idx[2] = tri.vertex_[2]; + + const vertex_t & A = grid_.mesh_[idx[0]]; + const vertex_t & B = grid_.mesh_[idx[1]]; + const vertex_t & C = grid_.mesh_[idx[2]]; + + // calculate partial derivatrives of wB and wC with respect to A, B, C: + double dw[2][6]; + { + const double & Ax = A.xy_[0]; + const double & Ay = A.xy_[1]; + const double & Bx = B.xy_[0]; + const double & By = B.xy_[1]; + const double & Cx = C.xy_[0]; + const double & Cy = C.xy_[1]; + const double & Px = P[0]; + const double & Py = P[1]; + + double bx = Bx - Ax; + double cx = Cx - Ax; + double px = Px - Ax; + + double by = By - Ay; + double cy = Cy - Ay; + double py = Py - Ay; + + double pycx_pxcy = py * cx - px * cy; + double pxby_pybx = px * by - py * bx; + double bycx_bxcy = by * cx - bx * cy; + double bycx_bxcy_2 = bycx_bxcy * bycx_bxcy; + + // dwB/dAx, dwB/dAy: + dw[0][0] = (Cy - Py) / bycx_bxcy - (Cy - By) / bycx_bxcy_2; + dw[0][1] = (Px - Cx) / bycx_bxcy - (Bx - Cx) / bycx_bxcy_2; + + // dwB/dBx, dwB/dBy: + dw[0][2] = cy * pycx_pxcy / bycx_bxcy_2; + dw[0][3] = -cx * pycx_pxcy / bycx_bxcy_2; + + // dwB/dCx, dwB/dCy: + dw[0][4] = py / bycx_bxcy - by * pycx_pxcy / bycx_bxcy_2; + dw[0][5] = bx * pycx_pxcy / bycx_bxcy_2 - px / bycx_bxcy; + + // dwC/dAx, dwC/dAy: + dw[1][0] = (Py - By) / bycx_bxcy - (Cy - By) * pxby_pybx / bycx_bxcy_2; + dw[1][1] = (Bx - Px) / bycx_bxcy - (Bx - Cx) * pxby_pybx / bycx_bxcy_2; + + // dwC/dBx, dwC/dBy: + dw[1][2] = cy * pxby_pybx / bycx_bxcy_2 - py / bycx_bxcy; + dw[1][3] = px / bycx_bxcy - cx * pxby_pybx / bycx_bxcy_2; + + // dwC/dCx, dwC/dCy: + dw[1][4] = -by * pxby_pybx / bycx_bxcy_2; + dw[1][5] = bx * pxby_pybx / bycx_bxcy_2; + } + + double bu = B.uv_[0] - A.uv_[0]; + double bv = B.uv_[1] - A.uv_[1]; + double cu = C.uv_[0] - A.uv_[0]; + double cv = C.uv_[1] - A.uv_[1]; + + // shortcut: + double * du = &(jac[0]); + double * dv = &(jac[6]); + + for (unsigned int i = 0; i < 6; i++) + { + du[i] = bu * dw[0][i] + cu * dw[1][i]; + dv[i] = bv * dw[0][i] + cv * dw[1][i]; + } + + return true; +} + + +//---------------------------------------------------------------- +// the_grid_transform_t::the_grid_transform_t +// +the_grid_transform_t::the_grid_transform_t() + : rows_(0) + , cols_(0) +{} + +//---------------------------------------------------------------- +// the_grid_transform_t::is_ready +// +bool +the_grid_transform_t::is_ready() const +{ + unsigned int n = rows_ * cols_; + return (n != 0) && (grid_.tri_.size() == 2 * n); +} + +//---------------------------------------------------------------- +// the_grid_transform_t::transform_inv +// +bool +the_grid_transform_t::transform_inv(const pnt2d_t & uv, pnt2d_t & xy) const +{ +#if 0 + const double & c = uv[0]; + int c0 = int(floor(c * double(cols_))); + int c1 = c0 + 1; + if (c0 < 0 || c0 > int(cols_)) return false; + if (c0 == int(cols_)) + { + c0--; + c1--; + } + + const double & r = uv[1]; + int r0 = int(floor(r * double(rows_))); + int r1 = r0 + 1; + if (r0 < 0 || r0 > int(rows_)) return false; + if (r0 == int(rows_)) + { + r0--; + r1--; + } + + double w[4]; + { + double wx[] = { + double(c1) - c * double(cols_), + c * double(cols_) - double(c0) + }; + + double wy[] = { + double(r1) - r * double(rows_), + r * double(rows_) - double(r0) + }; + + w[0] = wx[0] * wy[0]; + w[1] = wx[1] * wy[0]; + w[2] = wx[1] * wy[1]; + w[3] = wx[0] * wy[1]; + } + + xy[0] = + vertex(r0, c0).xy_[0] * w[0] + + vertex(r0, c1).xy_[0] * w[1] + + vertex(r1, c1).xy_[0] * w[2] + + vertex(r1, c0).xy_[0] * w[3]; + + xy[1] = + vertex(r0, c0).xy_[1] * w[0] + + vertex(r0, c1).xy_[1] * w[1] + + vertex(r1, c1).xy_[1] * w[2] + + vertex(r1, c0).xy_[1] * w[3]; + + return true; + +#else + + return the_base_triangle_transform_t::transform_inv(uv, xy); +#endif +} + +//---------------------------------------------------------------- +// the_grid_transform_t::setup +// +void +the_grid_transform_t::setup(unsigned int rows, + unsigned int cols, + const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const std::vector & xy) +{ + rows_ = rows; + cols_ = cols; + + tile_min_ = tile_min; + tile_ext_ = tile_max - tile_min; + + grid_.mesh_.resize((rows_ + 1) * (cols_ + 1)); + grid_.resize(rows_ * GRID_DENSITY, cols_ * GRID_DENSITY); + + // temporaries: + pnt2d_t uv; + + for (unsigned int row = 0; row <= rows_; row++) + { + uv[1] = double(row) / double(rows_); + for (unsigned int col = 0; col <= cols_; col++) + { + uv[0] = double(col) / double(cols_); + + vertex_t & vx = vertex(row, col); + vx.uv_ = uv; + vx.xy_ = xy[row * (cols_ + 1) + col]; + } + } + + setup_mesh(); +} + +//---------------------------------------------------------------- +// the_grid_transform_t::setup_mesh +// +void +the_grid_transform_t::setup_mesh() +{ + grid_.tri_.resize(rows_ * cols_ * 2); + + // temporaries: + unsigned int v_idx[4]; + unsigned int t_idx[2]; + + for (unsigned int row = 0; row < rows_; row++) + { + for (unsigned int col = 0; col < cols_; col++) + { + // vertex indices: + v_idx[0] = (cols_ + 1) * row + col; + v_idx[1] = v_idx[0] + (cols_ + 1); + v_idx[2] = v_idx[1] + 1; + v_idx[3] = v_idx[0] + 1; + + // triangle indices: + t_idx[0] = (cols_ * row + col) * 2; + t_idx[1] = t_idx[0] + 1; + + // setup triangle A: + triangle_t & tri_a = grid_.tri_[t_idx[0]]; + tri_a.vertex_[0] = v_idx[0]; + tri_a.vertex_[1] = v_idx[1]; + tri_a.vertex_[2] = v_idx[2]; + + // setup triangle B: + triangle_t & tri_b = grid_.tri_[t_idx[1]]; + tri_b.vertex_[0] = v_idx[0]; + tri_b.vertex_[1] = v_idx[2]; + tri_b.vertex_[2] = v_idx[3]; + } + } + + // initialize the acceleration grid: + grid_.rebuild(); +} + + +//---------------------------------------------------------------- +// the_mesh_transform_t::is_ready +// +bool +the_mesh_transform_t::is_ready() const +{ + bool ready = (grid_.tri_.size() > 1); + return ready; +} + +//---------------------------------------------------------------- +// the_mesh_transform_t::setup +// +bool +the_mesh_transform_t::setup(const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const std::vector & uv, + const std::vector & xy, + unsigned int accel_grid_rows, + unsigned int accel_grid_cols) +{ + const std::size_t num_pts = uv.size(); + + tile_min_ = tile_min; + tile_ext_ = tile_max - tile_min; + + grid_.mesh_.resize(num_pts); + grid_.resize(accel_grid_rows, accel_grid_cols); + + for (unsigned int i = 0; i < num_pts; i++) + { + vertex_t & vx = grid_.mesh_[i]; + vx.uv_ = uv[i]; + vx.xy_ = xy[i]; + } + + return setup_mesh(); +} + +//---------------------------------------------------------------- +// the_mesh_transform_t::insert_point +// +bool +the_mesh_transform_t::insert_point(const pnt2d_t & uv, const pnt2d_t & xy, const bool delay_setup) +{ + vertex_t vx; + vx.uv_ = uv; + vx.xy_ = xy; + + // Check for duplicate entries. + std::vector::iterator iter = grid_.mesh_.begin(); + for (; iter != grid_.mesh_.end(); ++iter) + { + if (vx.uv_ == (*iter).uv_ && vx.xy_ == (*iter).xy_) + break; + } + + // We can't allow duplicates, as it breaks the triangulation algorithm. + if (iter != grid_.mesh_.end()) + return false; + + grid_.mesh_.push_back(vx); + + // At times we want to add a bunch of points before performing + // Delauney. By waiting we can save cycles. + return (delay_setup) ? false : setup_mesh(); +} + +//---------------------------------------------------------------- +// the_mesh_transform_t::insert_point +// +bool +the_mesh_transform_t::insert_point(const pnt2d_t & uv) +{ + pnt2d_t xy; + if (!transform_inv(uv, xy)) + { + return false; + } + + return insert_point(uv, xy); +} + + +//---------------------------------------------------------------- +// CircumCircle +// +// Credit to Paul Bourke (pbourke@swin.edu.au) +// for the original Fortran 77 Program :)) +// +// Check out http://local.wasp.uwa.edu.au/~pbourke/papers/triangulate/index.html +// You can use this code however you like providing the above credits +// remain in tact. +// +// Return true if a point (xp, yp) is inside the circumcircle +// made up of the points (x1, y1), (x2, y2), (x3, y3) +// +// The circumcircle centre is returned in (xc,yc) and the radius r +// +// NOTE: A point on the edge is inside the circumcircle +// +static bool +CircumCircle(double xp, + double yp, + double x1, + double y1, + double x2, + double y2, + double x3, + double y3, + double * xc, + double * yc, + double * rsqr) +{ + double fabsy1y2 = fabs(y1 - y2); + double fabsy2y3 = fabs(y2 - y3); + + // Check for coincident points + if (fabsy1y2 < EPSILON && fabsy2y3 < EPSILON) + { + return (false); + } + + // temporaries: + double m1, m2, mx1, mx2, my1, my2; + + if (fabsy1y2 < EPSILON) + { + m2 = -(x3 - x2) / (y3 - y2); + mx2 = (x2 + x3) / 2.0; + my2 = (y2 + y3) / 2.0; + *xc = (x2 + x1) / 2.0; + *yc = m2 * (*xc - mx2) + my2; + } + else if (fabsy2y3 < EPSILON) + { + m1 = -(x2 - x1) / (y2 - y1); + mx1 = (x1 + x2) / 2.0; + my1 = (y1 + y2) / 2.0; + *xc = (x3 + x2) / 2.0; + *yc = m1 * (*xc - mx1) + my1; + } + else + { + m1 = -(x2 - x1) / (y2 - y1); + m2 = -(x3 - x2) / (y3 - y2); + mx1 = (x1 + x2) / 2.0; + mx2 = (x2 + x3) / 2.0; + my1 = (y1 + y2) / 2.0; + my2 = (y2 + y3) / 2.0; + *xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2); + + if (fabsy1y2 > fabsy2y3) + { + *yc = m1 * (*xc - mx1) + my1; + } + else + { + *yc = m2 * (*xc - mx2) + my2; + } + } + + double dx = x2 - *xc; + double dy = y2 - *yc; + *rsqr = dx * dx + dy * dy; + + dx = xp - *xc; + dy = yp - *yc; + + double drsqr = dx * dx + dy * dy; + return (drsqr <= *rsqr); +} + +//---------------------------------------------------------------- +// TVertex +// +typedef struct +{ + double x, y; + unsigned int idx; +} TVertex; + +//---------------------------------------------------------------- +// TTriangle +// +typedef struct +{ + int p1, p2, p3; +} TTriangle; + +//---------------------------------------------------------------- +// TEdge +// +typedef struct +{ + int p1, p2; +} TEdge; + +//---------------------------------------------------------------- +// vertex_sorter_t +// +static bool +TVertexCompare(const TVertex & a, const TVertex & b) +{ + // sort the vertices on the x-coordinate, in ascending order: + return a.x < b.x; +} + +//---------------------------------------------------------------- +// Triangulate +// +// Credit to Paul Bourke (pbourke@swin.edu.au) +// for the original Fortran 77 Program :)) +// +// Check out http://local.wasp.uwa.edu.au/~pbourke/papers/triangulate/index.html +// You can use this code however you like providing the above credits +// remain in tact. +// +// Takes as input a number of vertices in vertex array +// Passes back a list of triangular faces +// +// The triangles are arranged in a consistent clockwise order. +// The triangle array should be allocated to hold 2 * (nv + 3) - 2 triangles +// +// The vertex array must be big enough to hold 3 extra vertices +// (used internally for the supertiangle) +// +// The vertex array must be sorted in increasing x values +// +static int +Triangulate(int num_vertices, TVertex * vertex, TTriangle * tri, const int max_triangles, int & num_triangles) +{ + // Allocate memory for the completeness list, flag for each triangle + std::vector complete; + complete.assign(max_triangles, false); + + // Allocate memory for the edge list + int emax = 200; + std::vector edges(emax); + + // Set up the supertriangle + // + // This is a triangle which encompasses all the sample points. + // The supertriangle coordinates are added to the end of the + // vertex list. The supertriangle is the first triangle in + // the triangle list. + // + { + // Given a bounding box of all the vertices we can + // construct a bounding circle. Given the bounding circle we can + // construct a bounding triangle. + + // find the bounding box: + double xmin = vertex[0].x; + double ymin = vertex[0].y; + double xmax = xmin; + double ymax = ymin; + + for (int i = 1; i < num_vertices; i++) + { + if (vertex[i].x < xmin) + xmin = vertex[i].x; + if (vertex[i].x > xmax) + xmax = vertex[i].x; + if (vertex[i].y < ymin) + ymin = vertex[i].y; + if (vertex[i].y > ymax) + ymax = vertex[i].y; + } + + // find the (inflated) bounding circle: + double dx = xmax - xmin; + double dy = ymax - ymin; + + double xmid = (xmax + xmin) / 2.0; + double ymid = (ymax + ymin) / 2.0; + + double inflate = 1.1; + double r = (sqrt(dx * dx + dy * dy) / 2.0) * inflate; + + // find the bounding triangle with clockwise winding: + static const double sqrt3 = sqrt(3.0); + + vertex[num_vertices + 0].x = xmid - r * sqrt3; + vertex[num_vertices + 0].y = ymid - r; + vertex[num_vertices + 0].idx = num_vertices; + + vertex[num_vertices + 1].x = xmid; + vertex[num_vertices + 1].y = ymid + r * 2.0; + vertex[num_vertices + 1].idx = num_vertices + 1; + + vertex[num_vertices + 2].x = xmid + r * sqrt3; + vertex[num_vertices + 2].y = ymid - r; + vertex[num_vertices + 2].idx = num_vertices + 2; + + tri[0].p1 = num_vertices; + tri[0].p2 = num_vertices + 1; + tri[0].p3 = num_vertices + 2; + + complete[0] = false; + num_triangles = 1; + } + + // temporaries: + double xp = 0; + double yp = 0; + double x1 = 0; + double y1 = 0; + double x2 = 0; + double y2 = 0; + double x3 = 0; + double y3 = 0; + double xc = 0; + double yc = 0; + double r = 0; + + // Include each point one at a time into the existing mesh + int nedge = 0; + + for (int i = 0; i < num_vertices; i++) + { + xp = vertex[i].x; + yp = vertex[i].y; + nedge = 0; + + // Set up the edge buffer. + // If the point (xp,yp) lies inside the circumcircle then the + // three edges of that triangle are added to the edge buffer + // and that triangle is removed. + // + for (int j = 0; j < num_triangles; j++) + { + if (complete[j]) + { + continue; + } + + x1 = vertex[tri[j].p1].x; + y1 = vertex[tri[j].p1].y; + x2 = vertex[tri[j].p2].x; + y2 = vertex[tri[j].p2].y; + x3 = vertex[tri[j].p3].x; + y3 = vertex[tri[j].p3].y; + + const bool inside = CircumCircle(xp, yp, x1, y1, x2, y2, x3, y3, &xc, &yc, &r); + + if (xc < xp && ((xp - xc) * (xp - xc)) > r) + { + complete[j] = true; + } + + if (inside) + { + // Check that we haven't exceeded the edge list size + if (nedge + 3 >= emax) + { + emax += 100; + edges.resize(emax); + } + + edges[nedge + 0].p1 = tri[j].p1; + edges[nedge + 0].p2 = tri[j].p2; + edges[nedge + 1].p1 = tri[j].p2; + edges[nedge + 1].p2 = tri[j].p3; + edges[nedge + 2].p1 = tri[j].p3; + edges[nedge + 2].p2 = tri[j].p1; + + nedge += 3; + tri[j] = tri[num_triangles - 1]; + complete[j] = complete[num_triangles - 1]; + + num_triangles--; + j--; + } + } + + // Tag multiple edges + // + // NOTE: if all triangles are specified anticlockwise then all + // interior edges are opposite pointing in direction. + // + for (int j = 0; j < nedge - 1; j++) + { + for (int k = j + 1; k < nedge; k++) + { + if ((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1)) + { + edges[j].p1 = -1; + edges[j].p2 = -1; + edges[k].p1 = -1; + edges[k].p2 = -1; + } + + // Shouldn't need the following, see note above + if ((edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2)) + { + edges[j].p1 = -1; + edges[j].p2 = -1; + edges[k].p1 = -1; + edges[k].p2 = -1; + } + } + } + + // Form new triangles for the current point + // Skipping over any tagged edges. + // All edges are arranged in clockwise order. + // + for (int j = 0; j < nedge; j++) + { + if (edges[j].p1 < 0 || edges[j].p2 < 0) + { + continue; + } + + if (num_triangles >= max_triangles) + { + return 4; + } + + tri[num_triangles].p1 = edges[j].p1; + tri[num_triangles].p2 = edges[j].p2; + tri[num_triangles].p3 = i; + complete[num_triangles] = false; + num_triangles++; + } + } + + // Remove triangles with supertriangle vertices + // (triangles which have a vertex number greater than num_vertices) + // + for (int i = 0; i < num_triangles; i++) + { + if (tri[i].p1 >= num_vertices || tri[i].p2 >= num_vertices || tri[i].p3 >= num_vertices) + { + tri[i] = tri[num_triangles - 1]; + num_triangles--; + i--; + } + } + + return 0; +} + + +//---------------------------------------------------------------- +// the_mesh_transform_t::setup_mesh +// +bool +the_mesh_transform_t::setup_mesh() +{ + std::size_t num_vertices = grid_.mesh_.size(); + if (num_vertices == 0) + { + return false; + } + + // shortcut: + const vertex_t * verts = &(grid_.mesh_[0]); + + // added 3 extra points for the super-triangle (bounding triangle): + std::vector vertex_vec(num_vertices + 3); + for (std::size_t i = 0; i < num_vertices; i++) + { + TVertex & vertex = vertex_vec[i]; + + // triangulation happens in the uv-space: + vertex.x = verts[i].uv_[0]; + vertex.y = verts[i].uv_[1]; + + // keep track of the original vertex index: + vertex.idx = i; + } + + // sort the vertices on the x-coordinate, in ascending order: + std::sort(vertex_vec.begin(), vertex_vec.begin() + num_vertices, &TVertexCompare); + + // In 2-D the number of triangles is calculated as + // Nt = 2 * Nv - 2 - Nb, where Nv is the number of vertices + // and Nb is the number of boundary vertices (vertices on the convex hull). + // + // Since we don't know how many of the vertices are boundary vertices + // we allocate a triangle buffer larger than necessary + // also accounting for the supertriangle: + // + const int max_triangles = 2 * (num_vertices + 3) - 2; + std::vector triangle_vec(max_triangles); + + // shortcuts: + TVertex * vertices = &(vertex_vec[0]); + TTriangle * triangles = &(triangle_vec[0]); + + int num_triangles = 0; + int error = ::Triangulate(num_vertices, vertices, triangles, max_triangles, num_triangles); + if (error != 0) + { + // this shouldn't happen: + return false; + } + + if (num_triangles == 0) + { + return true; + } + + // setup the mesh triangles: + grid_.tri_.resize(num_triangles); + triangle_t * tris = &(grid_.tri_[0]); + + for (int i = 0; i < num_triangles; i++) + { + triangle_t & tri = tris[i]; + + // flip triangle winding to counterclockwise: + tri.vertex_[0] = (vertices[triangles[i].p1]).idx; + tri.vertex_[2] = (vertices[triangles[i].p2]).idx; + tri.vertex_[1] = (vertices[triangles[i].p3]).idx; + } + + // initialize the acceleration grid: + grid_.rebuild(); + + return true; +} diff --git a/src/IRMosaicRefinementCommon.cxx b/src/IRMosaicRefinementCommon.cxx new file mode 100644 index 0000000..7cb74b1 --- /dev/null +++ b/src/IRMosaicRefinementCommon.cxx @@ -0,0 +1,171 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : mosaic_refinement_common.cxx +// Author : Pavel A. Koshevoy +// Created : Mon Nov 3 20:26:25 MST 2008 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : Helper functions for automatic mosaic refinement. + +// local includes: +#include "itk/mosaic_refinement_common.hxx" + + +//---------------------------------------------------------------- +// regularize_displacements +// +void +regularize_displacements(// computed displacement vectors of the moving image + // grid transform control points, in mosaic space: + std::vector & xy_shift, + std::vector & mass, + + image_t::Pointer & dx, + image_t::Pointer & dy, + image_t::Pointer & db, + + // median filter radius: + const unsigned int & median_radius) +{ + // shortcuts: + image_t::RegionType::SizeType sz = dx->GetLargestPossibleRegion().GetSize(); + unsigned int mesh_cols = sz[0]; + unsigned int mesh_rows = sz[1]; + + // denoise + if (median_radius > 0) + { + dx = median(dx, median_radius); + dy = median(dy, median_radius); + // db = median(db, median_radius); + } + + // extend (fill in gaps): + typedef itk::ImageRegionConstIteratorWithIndex iter_t; + iter_t iter(dx, dx->GetLargestPossibleRegion()); + image_t::Pointer dx_blurred = cast(dx); + image_t::Pointer dy_blurred = cast(dy); + image_t::Pointer db_blurred = cast(db); + + for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter) + { + image_t::IndexType index = iter.GetIndex(); + if (!db->GetPixel(index)) + { + static const double max_w = 3.0; + double w = 0.0; + double px = 0.0; + double py = 0.0; + + // keep expanding the neighborhood until at least one + // successful sample is found: + + int max_x = std::max(int(index[0]), int(mesh_cols - 1 - index[0])); + int max_y = std::max(int(index[1]), int(mesh_rows - 1 - index[1])); + int max_r = std::min(1, std::max(max_x, max_y)); + for (int r = 1; r <= max_r && w < max_w; r++) + { + image_t::IndexType ix; + int x0 = index[0] - r; + int x1 = index[0] + r; + int y0 = index[1] - r; + int y1 = index[1] + r; + + int d = 2 * r + 1; + for (int o = 0; o < d; o++) + { + ix[0] = x0; + ix[1] = y0 + o + 1; + if (ix[0] >= 0 && ix[0] < int(mesh_cols) && + ix[1] >= 0 && ix[1] < int(mesh_rows) && + db->GetPixel(ix)) + { + px += dx->GetPixel(ix); + py += dy->GetPixel(ix); + w += 1.0; + } + + ix[0] = x1; + ix[1] = y0 + o; + if (ix[0] >= 0 && ix[0] < int(mesh_cols) && + ix[1] >= 0 && ix[1] < int(mesh_rows) && + db->GetPixel(ix)) + { + px += dx->GetPixel(ix); + py += dy->GetPixel(ix); + w += 1.0; + } + + ix[0] = x0 + o; + ix[1] = y0; + if (ix[0] >= 0 && ix[0] < int(mesh_cols) && + ix[1] >= 0 && ix[1] < int(mesh_rows) && + db->GetPixel(ix)) + { + px += dx->GetPixel(ix); + py += dy->GetPixel(ix); + w += 1.0; + } + + ix[0] = x0 + o + 1; + ix[1] = y1; + if (ix[0] >= 0 && ix[0] < int(mesh_cols) && + ix[1] >= 0 && ix[1] < int(mesh_rows) && + db->GetPixel(ix)) + { + px += dx->GetPixel(ix); + py += dy->GetPixel(ix); + w += 1.0; + } + } + } + + if (w != 0.0) + { + dx_blurred->SetPixel(index, px / w); + dy_blurred->SetPixel(index, py / w); + db_blurred->SetPixel(index, 1); + } + } + } + + // blur: + dx_blurred = smooth(dx_blurred, 1.0); + dy_blurred = smooth(dy_blurred, 1.0); + // db_blurred = smooth(db_blurred, 1.0); + + dx = dx_blurred; + dy = dy_blurred; + db = db_blurred; + + // update the mesh displacement field: + iter = iter_t(dx, dx->GetLargestPossibleRegion()); + for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter) + { + image_t::IndexType index = iter.GetIndex(); + unsigned int i = index[0] + index[1] * mesh_cols; + + xy_shift[i][0] = dx->GetPixel(index); + xy_shift[i][1] = dy->GetPixel(index); + mass[i] += db->GetPixel(index); + } +} diff --git a/src/itkRegularStepGradientDescentOptimizer2.cxx b/src/itkRegularStepGradientDescentOptimizer2.cxx new file mode 100644 index 0000000..4121609 --- /dev/null +++ b/src/itkRegularStepGradientDescentOptimizer2.cxx @@ -0,0 +1,335 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkRegularStepGradientDescentOptimizer2.cxx +// Author : Pavel A. Koshevoy, Tolga Tasdizen +// Created : 2005/11/11 14:54 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An enhanced version of the +// itk::RegularStepGradientDescentOptimizer +// fixing a bug with relaxation and adding support for +// step size increase (pick_up_pace), back tracking and +// keeping track of the best metric value +// and associated parameters. + +#ifndef _itkRegularStepGradientDescentOptimizer2_txx +#define _itkRegularStepGradientDescentOptimizer2_txx + +// local includes: +#include "itkRegularStepGradientDescentOptimizer2.h" + +// ITK includes: +#include +#include +#include + +namespace itk +{ +#if 0 +} +#endif + +//---------------------------------------------------------------- +// RegularStepGradientDescentOptimizer2::RegularStepGradientDescentOptimizer2 +// +RegularStepGradientDescentOptimizer2::RegularStepGradientDescentOptimizer2() + : m_Maximize(false) + , m_Value(0.0) + , m_PreviousValue(0.0) + , m_GradientMagnitudeTolerance(1e-4) + , m_MaximumStepLength(1.0) + , m_MinimumStepLength(1e-3) + , m_CurrentStepLength(0.0) + , m_RelaxationFactor(0.5) + , m_StopCondition(MaximumNumberOfIterations) + , m_NumberOfIterations(100) + , m_CurrentIteration(0) + , m_BestValue(0.0) + , m_BackTracking(false) + , m_PickUpPaceSteps(1000000) + , log_(cerr_log()) +{ + itkDebugMacro("Constructor"); + m_CostFunction = NULL; +} + +//---------------------------------------------------------------- +// RegularStepGradientDescentOptimizer2::StartOptimization +// +// Start the optimization +// +void +RegularStepGradientDescentOptimizer2::StartOptimization() +{ + itkDebugMacro("StartOptimization"); + const unsigned int spaceDimension = m_CostFunction->GetNumberOfParameters(); + + if (m_RelaxationFactor < 0.0) + { + itkExceptionMacro(<< "Relaxation factor must be positive. Current value is " << m_RelaxationFactor); + return; + } + + if (m_RelaxationFactor >= 1.0) + { + itkExceptionMacro(<< "Relaxation factor must less than 1.0. Current value is " << m_RelaxationFactor); + return; + } + + // Make sure the scales have been set properly + const unsigned int scalesSize = this->GetScales().size(); + if (scalesSize != spaceDimension) + { + itkExceptionMacro(<< "The size of Scales is " << scalesSize + << ", but the NumberOfParameters for the CostFunction is " << spaceDimension << "."); + } + + m_CurrentStepLength = m_MaximumStepLength; + m_CurrentIteration = 0; + + m_BestValue = this->m_Maximize ? -std::numeric_limits::max() : std::numeric_limits::max(); + m_Value = m_BestValue; + m_PreviousValue = m_BestValue; + + m_Gradient = DerivativeType(spaceDimension); + m_Gradient.Fill(0.0f); + + m_PreviousGradient = DerivativeType(spaceDimension); + m_PreviousGradient.Fill(0.0f); + + this->SetCurrentPosition(GetInitialPosition()); + this->ResumeOptimization(); +} + +//---------------------------------------------------------------- +// RegularStepGradientDescentOptimizer2::ResumeOptimization +// +// Resume the optimization +// +void +RegularStepGradientDescentOptimizer2::ResumeOptimization() +{ + // keep track of the number of sequential steps that + // resulted in function optimization: + unsigned int successful_steps = 0; + + itkDebugMacro("ResumeOptimization"); + this->InvokeEvent(StartEvent()); + + m_Stop = false; + while (!m_Stop) + { + m_PreviousGradient = m_Gradient; + m_PreviousValue = m_Value; + + ParametersType currentPosition = this->GetCurrentPosition(); + m_CostFunction->GetValueAndDerivative(currentPosition, m_Value, m_Gradient); + + if ((this->m_Maximize && m_Value < m_PreviousValue) || (!this->m_Maximize && m_Value > m_PreviousValue)) + { + // relax the step size: + m_CurrentStepLength *= m_RelaxationFactor; + + if (m_BackTracking) + { + // FIXME: + *log_ << m_Value << " vs " << m_BestValue + << " -- relaxing and backtracking, new step length: " << m_CurrentStepLength << std::endl; + + // backtrack to the previous position: + this->SetCurrentPosition(m_BestParams); + currentPosition = this->GetCurrentPosition(); + m_CostFunction->GetValueAndDerivative(currentPosition, m_Value, m_Gradient); + } + else + { + // FIXME: + *log_ << m_Value << " vs " << m_PreviousValue << " -- relaxing, new step length: " << m_CurrentStepLength + << std::endl; + } + + successful_steps = 0; + } + else + { + successful_steps++; + if (successful_steps % m_PickUpPaceSteps == 0) + { + // pick up the pace after N successful steps: + m_CurrentStepLength = std::min(m_MaximumStepLength, m_CurrentStepLength / m_RelaxationFactor); + + // FIXME: + *log_ << successful_steps << " successful steps -- increasing pace, new step length: " << m_CurrentStepLength + << std::endl; + } + } + + if (m_CurrentIteration == m_NumberOfIterations) + { + m_Stop = true; + m_StopCondition = MaximumNumberOfIterations; + this->StopOptimization(); + } + + if ((this->m_Maximize && m_Value > m_BestValue) || (!this->m_Maximize && m_Value < m_BestValue)) + { + m_BestValue = m_Value; + m_BestParams = currentPosition; + } + + if (!m_Stop) + { + this->AdvanceOneStep(); + m_CurrentIteration++; + } + } +} + +//---------------------------------------------------------------- +// RegularStepGradientDescentOptimizer2::StopOptimization +// +// Stop optimization +// +void +RegularStepGradientDescentOptimizer2::StopOptimization() +{ + itkDebugMacro("StopOptimization"); + m_Stop = true; + this->InvokeEvent(EndEvent()); +} + +//---------------------------------------------------------------- +// RegularStepGradientDescentOptimizer2::AdvanceOneStep +// +// Advance one Step following the gradient direction +// +void +RegularStepGradientDescentOptimizer2::AdvanceOneStep() +{ + itkDebugMacro("AdvanceOneStep"); + + const unsigned int spaceDimension = m_CostFunction->GetNumberOfParameters(); + + DerivativeType transformedGradient(spaceDimension); + DerivativeType previousTransformedGradient(spaceDimension); + const ScalesType & scales = this->GetScales(); + + for (unsigned int i = 0; i < spaceDimension; i++) + { + transformedGradient[i] = m_Gradient[i] / scales[i]; + previousTransformedGradient[i] = m_PreviousGradient[i] / scales[i]; + } + + double magnitudeSquared = 0.0; + for (unsigned int dim = 0; dim < spaceDimension; dim++) + { + const double weighted = transformedGradient[dim]; + magnitudeSquared += weighted * weighted; + } + + const double gradientMagnitude = vcl_sqrt(magnitudeSquared); + if (gradientMagnitude < m_GradientMagnitudeTolerance) + { + m_StopCondition = GradientMagnitudeTolerance; + this->StopOptimization(); + return; + } + + if (m_CurrentStepLength < m_MinimumStepLength) + { + m_StopCondition = StepTooSmall; + this->StopOptimization(); + return; + } + + const double direction = this->m_Maximize ? 1.0 : -1.0; + const double step_scale = direction * m_CurrentStepLength; + + // FIXME: + // *log_ << "gradientMagnitude: " << gradientMagnitude << std::endl; + + this->StepAlongGradient(step_scale, transformedGradient); + this->InvokeEvent(IterationEvent()); +} + +//---------------------------------------------------------------- +// RegularStepGradientDescentOptimizer2::StepAlongGradient +// +// Advance one Step following the gradient direction +// This method will be overridden in non-vector spaces +// +void +RegularStepGradientDescentOptimizer2::StepAlongGradient(double factor, const DerivativeType & transformedGradient) +{ + itkDebugMacro(<< "factor = " << factor << " transformedGradient = " << transformedGradient); + + const unsigned int spaceDimension = m_CostFunction->GetNumberOfParameters(); + + ParametersType newPosition(spaceDimension); + ParametersType currentPosition = this->GetCurrentPosition(); + + for (unsigned int j = 0; j < spaceDimension; j++) + { + newPosition[j] = currentPosition[j] + transformedGradient[j] * factor; + } + + itkDebugMacro(<< "new position = " << newPosition); + this->SetCurrentPosition(newPosition); +} + +//---------------------------------------------------------------- +// RegularStepGradientDescentOptimizer2::PrintSelf +// +void +RegularStepGradientDescentOptimizer2::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "MaximumStepLength: " << m_MaximumStepLength << std::endl; + os << indent << "MinimumStepLength: " << m_MinimumStepLength << std::endl; + os << indent << "RelaxationFactor: " << m_RelaxationFactor << std::endl; + os << indent << "GradientMagnitudeTolerance: " << m_GradientMagnitudeTolerance << std::endl; + os << indent << "NumberOfIterations: " << m_NumberOfIterations << std::endl; + os << indent << "CurrentIteration: " << m_CurrentIteration << std::endl; + os << indent << "Value: " << m_Value << std::endl; + os << indent << "Maximize: " << m_Maximize << std::endl; + + if (m_CostFunction) + { + os << indent << "CostFunction: " << &m_CostFunction << std::endl; + } + else + { + os << indent << "CostFunction: " << "(None)" << std::endl; + } + + os << indent << "CurrentStepLength: " << m_CurrentStepLength << std::endl; + os << indent << "StopCondition: " << m_StopCondition << std::endl; + os << indent << "Gradient: " << m_Gradient << std::endl; +} + +#if 0 +{ +#endif +} // end namespace itk + + +#endif // _itkRegularStepGradientDescentOptimizer2_txx