degate  0.1.2
TemplateMatching.cc
Go to the documentation of this file.
00001 /* -*-c++-*-
00002 
00003   This file is part of the IC reverse engineering tool degate.
00004 
00005   Copyright 2008, 2009, 2010 by Martin Schobert
00006 `
00007   Degate is free software: you can redistribute it and/or modify
00008   it under the terms of the GNU General Public License as published by
00009   the Free Software Foundation, either version 3 of the License, or
00010   any later version.
00011 
00012   Degate is distributed in the hope that it will be useful,
00013   but WITHOUT ANY WARRANTY; without even the implied warranty of
00014   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015   GNU General Public License for more details.
00016 
00017   You should have received a copy of the GNU General Public License
00018   along with degate. If not, see <http://www.gnu.org/licenses/>.
00019 
00020 */
00021 
00022 
00023 #include <TemplateMatching.h>
00024 #include <Image.h>
00025 #include <globals.h>
00026 #include <algorithm>
00027 #include <Statistics.h>
00028 #include <ImageHelper.h>
00029 #include <MedianFilter.h>
00030 #include <DegateHelper.h>
00031 
00032 #include <utility>
00033 #include <boost/foreach.hpp>
00034 #include <math.h>
00035 
00036 using namespace degate;
00037 
00038 //#define USE_MEDIAN_FILTER 2
00039 //#define USE_GAUSS_FILTER 5
00040 
00041 #if defined(USE_MEDIAN_FILTER) || defined(USE_GAUSS_FILTER)
00042 #define USE_FILTER
00043 #endif
00044 
00045 TemplateMatching::TemplateMatching() {
00046   threshold_hc = 0.40;
00047   threshold_detection = 0.70;
00048   max_step_size_search = 3;
00049   scale_down = 1;
00050 }
00051 
00052 TemplateMatching::~TemplateMatching() {
00053 }
00054 
00055 void TemplateMatching::precalc_sum_tables(TileImage_GS_BYTE_shptr img,
00056                                           TileImage_GS_DOUBLE_shptr summation_table_single,
00057                                           TileImage_GS_DOUBLE_shptr summation_table_squared) {
00058 
00059   unsigned int x, y;
00060 
00061   for(y = 0; y < img->get_height(); y++)
00062     for(x = 0; x < img->get_width(); x++) {
00063       summation_table_single->set_pixel(x, y, 0);
00064       summation_table_squared->set_pixel(x, y, 0);
00065     }
00066 
00067   for(y = 0; y < img->get_height(); y++)
00068     for(x = 0; x < img->get_width(); x++) {
00069 
00070       double f = img->get_pixel_as<gs_double_pixel_t>(x, y);
00071 
00072       double s_l = x > 0 ? summation_table_single->get_pixel(x - 1, y) : 0;
00073       double s_o = y > 0 ? summation_table_single->get_pixel(x, y - 1) : 0;
00074       double s_lo = x > 0 && y > 0 ? summation_table_single->get_pixel(x - 1, y - 1) : 0;
00075 
00076 
00077       double s2_l = x > 0 ? summation_table_squared->get_pixel(x - 1, y) : 0;
00078       double s2_o = y > 0 ? summation_table_squared->get_pixel(x, y - 1) : 0;
00079       double s2_lo = x > 0 && y > 0 ? summation_table_squared->get_pixel(x - 1, y - 1) : 0;
00080 
00081       double f1 = f + s_l + s_o - s_lo;
00082       double f2 = f*f + s2_l + s2_o - s2_lo;
00083 
00084       summation_table_single->set_pixel(x, y, f1);
00085       summation_table_squared->set_pixel(x, y, f2);
00086     }
00087 
00088 }
00089 
00090 
00091 
00092 void TemplateMatching::init(BoundingBox const& bounding_box, Project_shptr project) {
00093 
00094   assert(project != NULL);
00095 
00096   this->project = project;
00097   this->bounding_box = bounding_box;
00098 
00099   // limit bounding box
00100   if(this->bounding_box.get_max_x() + 1 > (int)project->get_width())
00101     this->bounding_box.set_max_x(LENGTH_TO_MAX(project->get_width()));
00102 
00103   if(this->bounding_box.get_max_y() + 1 > (int)project->get_height())
00104     this->bounding_box.set_max_y(LENGTH_TO_MAX(project->get_height()));
00105 
00106   ScalingManager_shptr sm = layer_matching->get_scaling_manager();
00107 
00108   debug(TM, "Prepare background.");
00109   prepare_background_images(sm, bounding_box, get_scaling_factor());
00110   debug(TM, "Prepare sum tabes.");
00111   prepare_sum_tables(gs_img_normal, gs_img_scaled);
00112 
00113   reset_progress();
00114 }
00115 
00116 
00117 BoundingBox TemplateMatching::get_scaled_bounding_box(BoundingBox const& bounding_box,
00118                                                       double scale_down) const {
00119   return BoundingBox(lrint(bounding_box.get_min_x() / scale_down),
00120                      lrint(bounding_box.get_max_x() / scale_down),
00121                      lrint(bounding_box.get_min_y() / scale_down),
00122                      lrint(bounding_box.get_max_y() / scale_down));
00123 }
00124 
00125 void TemplateMatching::prepare_background_images(ScalingManager_shptr sm,
00126                                                  BoundingBox const& bounding_box,
00127                                                  unsigned int scaling_factor) {
00128 
00129   // Get the normal background image and the scaled background image
00130   // These images are in RGBA format.
00131   const ScalingManager<BackgroundImage>::image_map_element i1 =
00132     sm->get_image(1);
00133   const ScalingManager<BackgroundImage>::image_map_element i2 =
00134     sm->get_image(scaling_factor);
00135 
00136   assert(i1.second != NULL);
00137   assert(i2.second != NULL);
00138   assert(i2.first  == get_scaling_factor());
00139 
00140   BackgroundImage_shptr img_normal = i1.second;
00141   BackgroundImage_shptr img_scaled = i2.second;
00142 
00143   // Create a greyscaled image for the normal
00144   // unscaled background image and the scaled version.
00145   BoundingBox scaled_bounding_box =
00146     get_scaled_bounding_box(bounding_box, scaling_factor);
00147 
00148 
00149   gs_img_normal = TileImage_GS_BYTE_shptr(new TileImage_GS_BYTE(bounding_box.get_width(),
00150                                                                 bounding_box.get_height()));
00151 
00152 #ifdef USE_FILTER
00153 
00154   TileImage_GS_BYTE_shptr tmp;
00155 
00156   tmp =  TileImage_GS_BYTE_shptr(new TileImage_GS_BYTE(bounding_box.get_width(),
00157                                                        bounding_box.get_height()));;
00158 
00159   extract_partial_image(tmp, img_normal, bounding_box);
00160 
00161   #ifdef USE_MEDIAN_FILTER
00162 
00163   median_filter(gs_img_normal, tmp, USE_MEDIAN_FILTER);
00164 
00165 
00166   #elif defined(USE_GAUSS_FILTER)
00167 
00168   int blur_kernel_size = USE_GAUSS_FILTER;
00169 
00170   std::shared_ptr<GaussianBlur>
00171     gaussian_blur_kernel(new GaussianBlur(blur_kernel_size, blur_kernel_size, 1.1));
00172 
00173   gaussian_blur_kernel->print();
00174 
00175   convolve(gs_img_normal, tmp, gaussian_blur_kernel);
00176 
00177   #endif
00178 
00179 #else
00180   extract_partial_image(gs_img_normal, img_normal, bounding_box);
00181 #endif
00182 
00183 
00184 
00185   if(scaling_factor == 1)
00186     gs_img_scaled = gs_img_normal;
00187   else {
00188 
00189     gs_img_scaled = TileImage_GS_BYTE_shptr(new TileImage_GS_BYTE(scaled_bounding_box.get_width(),
00190                                                                   scaled_bounding_box.get_height()));
00191 
00192 #ifdef USE_MEDIAN_FILTER
00193     tmp =  TileImage_GS_BYTE_shptr(new TileImage_GS_BYTE(bounding_box.get_width(),
00194                                                          bounding_box.get_height()));;
00195 
00196     extract_partial_image(tmp, img_scaled, scaled_bounding_box);
00197 
00198     median_filter(gs_img_scaled, tmp, USE_MEDIAN_FILTER);
00199 #else
00200 
00201     extract_partial_image(gs_img_scaled, img_scaled, scaled_bounding_box);
00202 #endif
00203   }
00204 
00205   //save_image("/tmp/xxx1.tif", gs_img_normal);
00206   //save_image("/tmp/xxx2.tif", gs_img_scaled);
00207 
00208 }
00209 
00210 void TemplateMatching::prepare_sum_tables(TileImage_GS_BYTE_shptr gs_img_normal,
00211                                           TileImage_GS_BYTE_shptr gs_img_scaled) {
00212 
00213   unsigned int
00214     w_n = gs_img_normal->get_width(),
00215     h_n = gs_img_normal->get_height(),
00216     w_s = gs_img_scaled->get_width(),
00217     h_s = gs_img_scaled->get_height();
00218 
00219   sum_table_single_normal = TileImage_GS_DOUBLE_shptr(new TileImage_GS_DOUBLE(w_n, h_n));
00220   sum_table_squared_normal = TileImage_GS_DOUBLE_shptr(new TileImage_GS_DOUBLE(w_n, h_n));
00221   sum_table_single_scaled = TileImage_GS_DOUBLE_shptr(new TileImage_GS_DOUBLE(w_s, h_s));
00222   sum_table_squared_scaled = TileImage_GS_DOUBLE_shptr(new TileImage_GS_DOUBLE(w_s, h_s));
00223 
00224   precalc_sum_tables(gs_img_normal, sum_table_single_normal, sum_table_squared_normal);
00225   precalc_sum_tables(gs_img_scaled, sum_table_single_scaled, sum_table_squared_scaled);
00226 }
00227 
00228 
00229 bool compare_template_size(const GateTemplate_shptr lhs, const GateTemplate_shptr rhs) {
00230   return lhs->get_width() * lhs->get_height() > rhs->get_width() * rhs->get_height();
00231 }
00232 
00233 bool compare_correlation(TemplateMatching::match_found const& lhs,
00234                          TemplateMatching::match_found const&  rhs) {
00235   return lhs.correlation > rhs.correlation;
00236 }
00237 
00238 void TemplateMatching::set_templates(std::list<GateTemplate_shptr> tmpl_set) {
00239   this->tmpl_set = tmpl_set;
00240   this->tmpl_set.sort(compare_template_size);
00241 }
00242 
00243 void TemplateMatching::set_orientations(std::list<Gate::ORIENTATION> tmpl_orientations) {
00244   this->tmpl_orientations = tmpl_orientations;
00245   for(std::list<Gate::ORIENTATION>::const_iterator iter = tmpl_orientations.begin();
00246       iter != tmpl_orientations.end(); ++iter) {
00247     debug(TM, "Template orientation to match: %d", *iter);
00248   }
00249 }
00250 
00251 void TemplateMatching::run() {
00252 
00253 
00254   debug(TM, "run template matching");
00255   std::list<match_found> matches;
00256 
00257   stats.reset();
00258   set_progress_step_size(1.0/(tmpl_set.size() * tmpl_orientations.size()));
00259 
00260   /*
00261   BOOST_FOREACH(GateTemplate_shptr tmpl, tmpl_set) {
00262     BOOST_FOREACH(Gate::ORIENTATION orientation, tmpl_orientations) {
00263       add_task( run2, parameter)
00264     }
00265   }
00266 
00267   wait;
00268   */
00269   BOOST_FOREACH(GateTemplate_shptr tmpl, tmpl_set) {
00270 
00271     BOOST_FOREACH(Gate::ORIENTATION orientation, tmpl_orientations) {
00272 
00273       boost::format f("Check cell \"%1%\"");
00274       f % tmpl->get_name();
00275       set_log_message(f.str());
00276 
00277       prepared_template prep_tmpl_img = prepare_template(tmpl, orientation);
00278       //match_single_template(prep_tmpl_img, t_hc, t_det);
00279 
00280 
00281       std::list<match_found> m = match_single_template(prep_tmpl_img,
00282                                                        threshold_hc,
00283                                                        threshold_detection);
00284 
00285 
00286       matches.insert(matches.end(), m.begin(), m.end());
00287 
00288       progress_step_done();
00289       if(is_canceled()) {
00290         reset_progress();
00291         return;
00292       }
00293     }
00294 
00295   }
00296 
00297 
00298   matches.sort(compare_correlation);
00299 
00300   BOOST_FOREACH(match_found const& m, matches) {
00301     std::cout << "Try to insert gate of type " << m.tmpl->get_name() << " with corr="
00302               << m.correlation << " at " << m.x << "," << m.y << std::endl;
00303     if(add_gate(m.x, m.y, m.tmpl, m.orientation, m.correlation, m.t_hc))
00304       std::cout << "\tInserted gate of type " << m.tmpl->get_name() << std::endl;
00305   }
00306 
00307   reset_progress();
00308 }
00309 
00310 
00311 double TemplateMatching::subtract_mean(TempImage_GS_BYTE_shptr img,
00312                                        TempImage_GS_DOUBLE_shptr zero_mean_img) const {
00313 
00314   double mean = average(img);
00315 
00316   double sum_over_zero_mean_img = 0;
00317   unsigned int x, y;
00318 
00319   for(y = 0; y < img->get_height(); y++)
00320     for(x = 0; x < img->get_width(); x++) {
00321       double tmp = img->get_pixel_as<gs_double_pixel_t>(x, y) - mean;
00322       zero_mean_img->set_pixel(x, y, tmp);
00323       sum_over_zero_mean_img += tmp * tmp;
00324     }
00325 
00326 
00327   return sum_over_zero_mean_img;
00328 }
00329 
00330 TemplateMatching::prepared_template TemplateMatching::prepare_template(GateTemplate_shptr tmpl,
00331                                                                        Gate::ORIENTATION orientation) {
00332 
00333   prepared_template prep;
00334 
00335   assert(layer_matching->get_layer_type() != Layer::UNDEFINED);
00336   assert(tmpl->has_image(layer_matching->get_layer_type()));
00337 
00338   prep.gate_template = tmpl;
00339   prep.orientation = orientation;
00340 
00341   // get image from template
00342   GateTemplateImage_shptr tmpl_img_orig = tmpl->get_image(layer_matching->get_layer_type());
00343 
00344   unsigned int
00345     w = tmpl_img_orig->get_width(),
00346     h = tmpl_img_orig->get_height();
00347 
00348   // get image according to orientation
00349   GateTemplateImage_shptr tmpl_img(new GateTemplateImage(w, h));
00350 
00351   //#ifdef USE_MEDIAN_FILERX
00352   //  median_filter(tmpl_img, tmpl_img_orig, 3);
00353   //#else
00354   copy_image(tmpl_img, tmpl_img_orig);
00355   //#endif
00356 
00357   //save_image("/tmp/xxx3.tif", tmpl_img);
00358 
00359   switch(orientation) {
00360   case Gate::ORIENTATION_FLIPPED_UP_DOWN:
00361     flip_up_down(tmpl_img);
00362     break;
00363   case Gate::ORIENTATION_FLIPPED_LEFT_RIGHT:
00364     flip_left_right(tmpl_img);
00365     break;
00366   case Gate::ORIENTATION_FLIPPED_BOTH:
00367     flip_both(tmpl_img);
00368     break;
00369   case Gate::ORIENTATION_NORMAL:
00370     break;
00371   case Gate::ORIENTATION_UNDEFINED:
00372     assert(1 == 0); // it is already checked
00373   }
00374 
00375   // create a scaled version of the template image
00376 
00377   unsigned int
00378     scaled_tmpl_width = (double)w / get_scaling_factor(),
00379     scaled_tmpl_height = (double)h / get_scaling_factor();
00380 
00381   prep.tmpl_img_normal = TempImage_GS_BYTE_shptr(new TempImage_GS_BYTE(w, h));
00382   copy_image(prep.tmpl_img_normal, tmpl_img);
00383 
00384   prep.tmpl_img_scaled = TempImage_GS_BYTE_shptr(new TempImage_GS_BYTE(scaled_tmpl_width,
00385                                                                        scaled_tmpl_height));
00386 
00387   scale_down_by_power_of_2(prep.tmpl_img_scaled, tmpl_img);
00388 
00389 
00390   // create zero-mean templates
00391   prep.zero_mean_template_normal = TempImage_GS_DOUBLE_shptr(new TempImage_GS_DOUBLE(w, h));
00392   prep.zero_mean_template_scaled = TempImage_GS_DOUBLE_shptr(new TempImage_GS_DOUBLE(scaled_tmpl_width,
00393                                                                                      scaled_tmpl_height));
00394 
00395 
00396   // subtract mean
00397 
00398   prep.sum_over_zero_mean_template_normal = subtract_mean(prep.tmpl_img_normal,
00399                                                           prep.zero_mean_template_normal);
00400   prep.sum_over_zero_mean_template_scaled = subtract_mean(prep.tmpl_img_scaled,
00401                                                           prep.zero_mean_template_scaled);
00402 
00403 
00404   assert(prep.sum_over_zero_mean_template_normal > 0);
00405   assert(prep.sum_over_zero_mean_template_scaled > 0);
00406 
00407   return prep;
00408 }
00409 
00410 
00411 
00412 void TemplateMatching::adjust_step_size(struct search_state & state, double corr_val) const {
00413   if(corr_val > 0) {
00414     state.step_size_search = std::max(1.0, rint((1.0 - (double)get_max_step_size()) *
00415                                                 corr_val + get_max_step_size()));
00416   }
00417 
00418   else state.step_size_search = get_max_step_size();
00419 }
00420 
00421 TemplateMatching::match_found
00422 TemplateMatching::keep_gate_match(unsigned int x, unsigned int y,
00423                                   struct prepared_template & tmpl,
00424                                   double corr_val, double threshold_hc) const {
00425   match_found hit;
00426   hit.x = x;
00427   hit.y = y;
00428   hit.tmpl = tmpl.gate_template;
00429   hit.correlation = corr_val;
00430   hit.t_hc = threshold_hc;
00431   hit.orientation = tmpl.orientation;
00432   return hit;
00433 }
00434 
00435 bool TemplateMatching::add_gate(unsigned int x, unsigned int y,
00436                                 GateTemplate_shptr tmpl,
00437                                 Gate::ORIENTATION orientation,
00438                                 double corr_val, double threshold_hc) {
00439 
00440   if(!layer_insert->exists_type_in_region<Gate>(x, x + tmpl->get_width(),
00441                                                 y, y + tmpl->get_height())) {
00442 
00443     Gate_shptr gate(new Gate(x, x + tmpl->get_width(),
00444                              y, y + tmpl->get_height(),
00445                              orientation));
00446 
00447     char dsc[100];
00448     snprintf(dsc, sizeof(dsc), "matched with corr=%.2f t_hc=%.2f", corr_val, threshold_hc);
00449     gate->set_description(dsc);
00450 
00451     gate->set_gate_template(tmpl);
00452 
00453     LogicModel_shptr lmodel = project->get_logic_model();
00454     lmodel->add_object(layer_insert, gate);
00455     lmodel->update_ports(gate);
00456 
00457     stats.hits++;
00458     return true;
00459   }
00460   return false;
00461 }
00462 
00463 std::list<TemplateMatching::match_found>
00464 TemplateMatching::match_single_template(struct prepared_template & tmpl,
00465                                         double threshold_hc, double threshold_detection) {
00466 
00467   debug(TM, "match_single_template(): start iterating over background image");
00468   search_state state;
00469   memset(&state, 0, sizeof(search_state));
00470   state.x = 1;
00471   state.y = 1;
00472   state.step_size_search = get_max_step_size();
00473   state.search_area = bounding_box;
00474   std::list<match_found> matches;
00475 
00476   double max_corr_for_search = -1;
00477 
00478   do { // works on unscaled, but cropped image
00479 
00480     double corr_val = calc_single_xcorr(gs_img_scaled,
00481                                         sum_table_single_scaled,
00482                                         sum_table_squared_scaled,
00483                                         tmpl.zero_mean_template_scaled,
00484                                         tmpl.sum_over_zero_mean_template_scaled,
00485                                         lrint((double)state.x / get_scaling_factor()),
00486                                         lrint((double)state.y / get_scaling_factor()));
00487 
00488     /*
00489     debug(TM, "%d,%d  == %d,%d  -> %f", state.x, state.y,
00490           lrint((double)state.x / get_scaling_factor()),
00491           lrint((double)state.y / get_scaling_factor()),
00492           corr_val);
00493     */
00494     if(corr_val > max_corr_for_search) max_corr_for_search = corr_val;
00495 
00496 
00497     adjust_step_size(state, corr_val);
00498 
00499     if(corr_val >= threshold_hc) {
00500       //debug(TM, "start hill climbing at(%d,%d), corr=%f", state.x, state.y, corr_val);
00501       unsigned int max_corr_x, max_corr_y;
00502       double curr_max_val;
00503       hill_climbing(state.x, state.y, corr_val,
00504                     &max_corr_x, &max_corr_y, &curr_max_val,
00505                     gs_img_normal, tmpl.zero_mean_template_normal,
00506                     tmpl.sum_over_zero_mean_template_normal);
00507 
00508       //debug(TM, "hill climbing returned for (%d,%d) corr=%f", max_corr_x, max_corr_y, curr_max_val);
00509       if(curr_max_val >= threshold_detection) {
00510         matches.push_back(keep_gate_match(max_corr_x + bounding_box.get_min_x(),
00511                                           max_corr_y + bounding_box.get_min_y(),
00512                                           tmpl, curr_max_val, threshold_hc));
00513       }
00514 
00515     }
00516 
00517   } while(get_next_pos(&state, tmpl) && !is_canceled());
00518 
00519   std::cout << "The maximum correlation value for the current template and orientation is " << max_corr_for_search << std::endl;
00520 
00521   return matches;
00522 }
00523 
00524 
00525 void TemplateMatching::hill_climbing(unsigned int start_x, unsigned int start_y, double xcorr_val,
00526                                      unsigned int * max_corr_x_out,
00527                                      unsigned int * max_corr_y_out,
00528                                      double * max_xcorr_out,
00529                                      const TileImage_GS_BYTE_shptr master,
00530                                      const TempImage_GS_DOUBLE_shptr zero_mean_template,
00531                                      double sum_over_zero_mean_template) const {
00532 
00533   unsigned int max_corr_x = start_x;
00534   unsigned int max_corr_y = start_y;
00535 
00536   double max_corr = xcorr_val;
00537   bool running = true;
00538 
00539   //std::list<std::pair<unsigned int, unsigned int> > positions;
00540 
00541   const unsigned int radius = get_max_step_size();
00542   const unsigned int size = (2 * radius + 1) * (2 * radius + 1);
00543 
00544   unsigned int positions_x[size], positions_y[size];
00545 
00546   while(running) {
00547     running = false;
00548 
00549     // first generate positions
00550 
00551   unsigned int
00552     from_x = max_corr_x >= radius ? max_corr_x - radius : 0,
00553     from_y = max_corr_y >= radius ? max_corr_y - radius : 0,
00554     to_x = max_corr_x + radius < master->get_width() ? max_corr_x + radius : master->get_width(),
00555     to_y = max_corr_y + radius < master->get_height() ? max_corr_y + radius : master->get_height();
00556 
00557     unsigned int i = 0;
00558     for(unsigned int _y = from_y; _y < to_y; _y++)
00559       for(unsigned int _x = from_x; _x < to_x; _x++) {
00560 
00561         if(max_corr_x != _x && max_corr_y != _y) {
00562           //positions.push_back(std::pair<unsigned int, unsigned int>(_x, _y));
00563           positions_x[i] = _x;
00564           positions_y[i] = _y;
00565           i++;
00566         }
00567       }
00568 
00569     // check for position with highest correlation value
00570     for(unsigned int i2 = 0; i2 < i; i2++) {
00571 
00572       unsigned int x = positions_x[i2], y = positions_y[i2];
00573 
00574       //debug(TM, "hill climbing step at (%d,%d)", x, y);
00575 
00576       double curr_corr_val = calc_single_xcorr(master,
00577                                                sum_table_single_normal,
00578                                                sum_table_squared_normal,
00579                                                zero_mean_template,
00580                                                sum_over_zero_mean_template,
00581                                                x, y);
00582 
00583       if(curr_corr_val > max_corr) {
00584         max_corr_x = x;
00585         max_corr_y = y;
00586         max_corr = curr_corr_val;
00587         running = true;
00588         //positions.clear();
00589       }
00590     }
00591   }
00592 
00593   *max_corr_x_out = max_corr_x;
00594   *max_corr_y_out = max_corr_y;
00595   *max_xcorr_out = max_corr;
00596 
00597 }
00598 
00599 
00600 double TemplateMatching::calc_single_xcorr(const TileImage_GS_BYTE_shptr master,
00601                                            const TileImage_GS_DOUBLE_shptr summation_table_single,
00602                                            const TileImage_GS_DOUBLE_shptr summation_table_squared,
00603                                            const TempImage_GS_DOUBLE_shptr zero_mean_template,
00604                                            double sum_over_zero_mean_template,
00605                                            unsigned int local_x,
00606                                            unsigned int local_y) const {
00607 
00608   double template_size = zero_mean_template->get_width() * zero_mean_template->get_height();
00609   assert(zero_mean_template->get_width() > 0 && zero_mean_template->get_height() > 0);
00610 
00611   unsigned int
00612     x_plus_w = local_x + zero_mean_template->get_width() -1,
00613     y_plus_h = local_y + zero_mean_template->get_height() -1,
00614     lxm1 = local_x - 1, // can wrap, it's checked later
00615     lym1 = local_y - 1;
00616 
00617   // calculate denominator
00618   double
00619     f1 = summation_table_single->get_pixel(x_plus_w, y_plus_h),
00620     f2 = summation_table_squared->get_pixel(x_plus_w, y_plus_h);
00621 
00622   if(local_x > 0) {
00623     f1 -= summation_table_single->get_pixel(lxm1, y_plus_h);
00624     f2 -= summation_table_squared->get_pixel(lxm1, y_plus_h);
00625   }
00626   if(local_y > 0) {
00627     f1 -= summation_table_single->get_pixel(x_plus_w, lym1);
00628     f2 -= summation_table_squared->get_pixel(x_plus_w, lym1);
00629   }
00630   if(local_x > 0 && local_y > 0) {
00631     f1 += summation_table_single->get_pixel(lxm1, lym1);
00632     f2 += summation_table_squared->get_pixel(lxm1, lym1);
00633   }
00634 
00635   double denominator = sqrt((f2 - f1*f1/template_size) * sum_over_zero_mean_template);
00636 
00637   // calculate nummerator
00638   if(std::isinf(denominator) || std::isnan(denominator) || denominator == 0) {
00639     debug(TM,
00640           "ERROR: The denominator is not a valid number: f1=%f f2=%f template_size=%f sum=%f "
00641           "local_x=%d local_y=%d x_plus_w=%d y_plus_h=%d lxm1=%d lym1=%d",
00642           f1, f2, template_size, sum_over_zero_mean_template,
00643           local_x, local_y, x_plus_w, y_plus_h, lxm1, lym1);
00644     return -1.0;
00645   }
00646 
00647   unsigned int _x, _y;
00648   double nummerator = 0;
00649 
00650   for(_y = 0; _y < zero_mean_template->get_height(); _y ++) {
00651     for(_x = 0; _x < zero_mean_template->get_width(); _x ++) {
00652       double f_xy = master->get_pixel(_x + local_x, _y + local_y);
00653       double t_xy = zero_mean_template->get_pixel(_x, _y);
00654       nummerator += f_xy * t_xy;
00655     }
00656   }
00657 
00658   double q = nummerator/denominator;
00659 
00660   //debug(TM, "x=%d,y=%d a=%f, nummerator=%f, denominator=%f, q=%f", local_x, local_y, a, nummerator, denominator, q);
00661 
00662   if(!(q >= -1.000001 && q <= 1.000001)) {
00663     debug(TM, "nummerator = %f / denominator = %f", nummerator, denominator);
00664   }
00665   assert(q >= -1.1 && q <= 1.1);
00666   return q;
00667 }
00668 
00669 bool TemplateMatchingNormal::get_next_pos(struct search_state * state,
00670                                           struct prepared_template const& tmpl) const {
00671 
00672   unsigned int
00673     tmpl_w = tmpl.tmpl_img_normal->get_width(),
00674     tmpl_h = tmpl.tmpl_img_normal->get_height();
00675 
00676   if(state->search_area.get_width() < tmpl_w ||
00677      state->search_area.get_height() < tmpl_h) return false;
00678 
00679   unsigned int step = state->step_size_search;
00680 
00681   bool there_was_a_gate = false;
00682 
00683   do {
00684     if(state->x + step < state->search_area.get_width() - tmpl_w)
00685       state->x += step;
00686     else {
00687       state->x = 1;
00688       if(state->y + step < state->search_area.get_height() - tmpl_h)
00689         state->y += step;
00690       else return false;
00691     }
00692 
00693     unsigned int dist_x =
00694       layer_insert->get_distance_to_gate_boundary(state->x + state->search_area.get_min_x(),
00695                                                   state->y + state->search_area.get_min_y(),
00696                                                   true, tmpl_w, tmpl_h);
00697 
00698     if(dist_x > 0) {
00699       /*
00700       debug(TM, "In the window starting at %d,%d there is already a gate. Skipping %d horizontal pixels",
00701             state->x + state->search_area.get_min_x(),
00702             state->y + state->search_area.get_min_y(),
00703             dist_x);
00704       */
00705       state->x += dist_x;
00706       there_was_a_gate = true;
00707       step = 1;
00708     }
00709     else there_was_a_gate = false;
00710 
00711   }while(there_was_a_gate);
00712 
00713   return true;
00714 }
00715 
00716 
00717 bool TemplateMatchingAlongGrid::initialize_state_struct(struct search_state * state,
00718                                                         int offs_min,
00719                                                         int offs_max,
00720                                                         bool is_horizontal_grid) const {
00721   if(state->grid == NULL) {
00722     const RegularGrid_shptr rg = is_horizontal_grid ?
00723       project->get_regular_horizontal_grid() : project->get_regular_vertical_grid();
00724     const IrregularGrid_shptr ig = is_horizontal_grid ?
00725       project->get_irregular_horizontal_grid() : project->get_irregular_vertical_grid();
00726 
00727     if(rg->is_enabled()) state->grid = rg;
00728     else if(ig->is_enabled()) state->grid = ig;
00729 
00730     if(state->grid != NULL) {
00731 
00732       debug(TM, "check grid from %d to %d", offs_min, offs_max);
00733 
00734       // iterate over grid and find first and last offset
00735       state->iter_begin = state->grid->begin();
00736       state->iter_last = state->grid->begin();
00737       state->iter_end = state->grid->end();
00738 
00739       for(Grid::grid_iter iter = state->grid->begin();
00740           iter != state->grid->end(); ++iter) {
00741         debug(TM, "\tcheck %d", *iter);
00742         if(*iter < offs_min) state->iter_begin = iter;
00743         if(*iter < offs_max) state->iter_last = iter;
00744       }
00745 
00746       if(*(state->iter_begin) < offs_min &&
00747          state->iter_begin != state->grid->end()) state->iter_begin++;
00748 
00749       //if(state->iter_last != state->grid->end()) state->iter_last++;
00750 
00751       if(state->iter_begin != state->grid->end())
00752         debug(TM, "first grid offset %d", *(state->iter_begin) );
00753       //      if(state->iter_last != state->grid->end())
00754       debug(TM, "last  grid offset %d", *(state->iter_last));
00755 
00756       if(state->iter_begin == state->grid->end() ||
00757          state->iter_last == state->grid->end()) {
00758         debug(TM, "failed");
00759         return false;
00760       }
00761 
00762       state->iter = state->iter_begin;
00763     }
00764     else {
00765       debug(TM, "There is no grid enabled.");
00766       return false;
00767     }
00768   }
00769   return true;
00770 }
00771 
00772 bool TemplateMatchingInRows::get_next_pos(struct search_state * state,
00773                                           struct prepared_template const& tmpl) const {
00774 
00775   // check if the search area is larger then the template
00776   unsigned int
00777     tmpl_w = tmpl.tmpl_img_normal->get_width(),
00778     tmpl_h = tmpl.tmpl_img_normal->get_height();
00779 
00780   if(state->search_area.get_width() < tmpl_w ||
00781      state->search_area.get_height() < tmpl_h) return false;
00782 
00783   unsigned int step = state->step_size_search;
00784 
00785   // get grid and check if we are working on regular or irregular grid
00786   if(state->grid == NULL &&
00787      initialize_state_struct(state,
00788                              state->search_area.get_min_y(),
00789                              state->search_area.get_max_y() - tmpl_h,
00790                              false) == false) {
00791     debug(TM, "Can't initialize search structure.");
00792     return false;
00793   }
00794 
00795 
00796   if(state->x == 1 && state->y == 1) { // start condition
00797     state->y = *(state->iter) - state->search_area.get_min_y();
00798   }
00799 
00800   bool there_was_a_gate = false;
00801 
00802   do {
00803 
00804     if(state->x + step < state->search_area.get_width() - tmpl_w)
00805       state->x += step;
00806     else {
00807       ++(state->iter);
00808 
00809       if(state->iter == state->iter_end ||
00810          *(state->iter) > *(state->iter_last)) {
00811         return false;
00812       }
00813 
00814       state->x = 1;
00815       state->y = *(state->iter) - state->search_area.get_min_y();
00816 
00817     }
00818 
00819     unsigned int dist_x = layer_insert->get_distance_to_gate_boundary(state->x + state->search_area.get_min_x(),
00820                                                                       state->y + state->search_area.get_min_y(),
00821                                                                       true, tmpl_w, tmpl_h);
00822 
00823     if(dist_x > 0) {
00824       debug(TM, "In the window starting at %d,%d there is already a gate. Skipping %d horizontal pixels",
00825             state->x + state->search_area.get_min_x(),
00826             state->y + state->search_area.get_min_y(),
00827             dist_x);
00828       state->x += dist_x;
00829       there_was_a_gate = true;
00830       step = 1;
00831     }
00832     else there_was_a_gate = false;
00833 
00834   }while(there_was_a_gate);
00835 
00836 
00837   return true;
00838 }
00839 
00840 
00841 bool TemplateMatchingInCols::get_next_pos(struct search_state * state,
00842                                           struct prepared_template const& tmpl) const {
00843 
00844   // check if the search area is larger then the template
00845   unsigned int
00846     tmpl_w = tmpl.tmpl_img_normal->get_width(),
00847     tmpl_h = tmpl.tmpl_img_normal->get_height();
00848 
00849   if(state->search_area.get_width() < tmpl_w ||
00850      state->search_area.get_height() < tmpl_h) return false;
00851 
00852   unsigned int step = state->step_size_search;
00853 
00854   // get grid and check if we are working on regular or irregular grid
00855   if(state->grid == NULL &&
00856      initialize_state_struct(state,
00857                              state->search_area.get_min_x(),
00858                              state->search_area.get_max_x() - tmpl_w,
00859                              true) == false)
00860     return false;
00861 
00862 
00863   if(state->x == 1 && state->y == 1) { // start condition
00864     state->x = *(state->iter) - state->search_area.get_min_x();
00865   }
00866 
00867   bool there_was_a_gate = false;
00868 
00869   do {
00870     if(state->y + step < state->search_area.get_height() - tmpl_h)
00871       state->y += step;
00872     else {
00873       ++state->iter;
00874 
00875       if(state->iter == state->iter_end ||
00876          *(state->iter) > *(state->iter_last)) return false;
00877 
00878       state->x = *(state->iter) - state->search_area.get_min_x();
00879       state->y = 1;
00880 
00881     }
00882 
00883     unsigned int dist_y = layer_insert->get_distance_to_gate_boundary(state->x + state->search_area.get_min_x(),
00884                                                                       state->y + state->search_area.get_min_y(),
00885                                                                       false, tmpl_w, tmpl_h);
00886 
00887     if(dist_y > 0) {
00888       debug(TM, "In the window starting at %d,%d there is already a gate. Skipping %d vertical pixels",
00889             state->x + state->search_area.get_min_x(),
00890             state->y + state->search_area.get_min_y(),
00891             dist_y);
00892       state->y += dist_y;
00893       there_was_a_gate = true;
00894       step = 1;
00895     }
00896     else there_was_a_gate = false;
00897 
00898   }while(there_was_a_gate);
00899 
00900   return true;
00901 }
00902