Home | History | Annotate | Download | only in kaze 
      1 //=============================================================================
      2 //
      3 // fed.cpp
      4 // Authors: Pablo F. Alcantarilla (1), Jesus Nuevo (2)
      5 // Institutions: Georgia Institute of Technology (1)
      6 //               TrueVision Solutions (2)
      7 // Date: 15/09/2013
      8 // Email: pablofdezalc (at) gmail.com
      9 //
     10 // AKAZE Features Copyright 2013, Pablo F. Alcantarilla, Jesus Nuevo
     11 // All Rights Reserved
     12 // See LICENSE for the license information
     13 //=============================================================================
     14 
     15 /**
     16  * @file fed.cpp
     17  * @brief Functions for performing Fast Explicit Diffusion and building the
     18  * nonlinear scale space
     19  * @date Sep 15, 2013
     20  * @author Pablo F. Alcantarilla, Jesus Nuevo
     21  * @note This code is derived from FED/FJ library from Grewenig et al.,
     22  * The FED/FJ library allows solving more advanced problems
     23  * Please look at the following papers for more information about FED:
     24  * [1] S. Grewenig, J. Weickert, C. Schroers, A. Bruhn. Cyclic Schemes for
     25  * PDE-Based Image Analysis. Technical Report No. 327, Department of Mathematics,
     26  * Saarland University, Saarbrcken, Germany, March 2013
     27  * [2] S. Grewenig, J. Weickert, A. Bruhn. From box filtering to fast explicit diffusion.
     28  * DAGM, 2010
     29  *
     30 */
     31 #include "../precomp.hpp"
     32 #include "fed.h"
     33 
     34 using namespace std;
     35 
     36 //*************************************************************************************
     37 //*************************************************************************************
     38 
     39 /**
     40  * @brief This function allocates an array of the least number of time steps such
     41  * that a certain stopping time for the whole process can be obtained and fills
     42  * it with the respective FED time step sizes for one cycle
     43  * The function returns the number of time steps per cycle or 0 on failure
     44  * @param T Desired process stopping time
     45  * @param M Desired number of cycles
     46  * @param tau_max Stability limit for the explicit scheme
     47  * @param reordering Reordering flag
     48  * @param tau The vector with the dynamic step sizes
     49  */
     50 int fed_tau_by_process_time(const float& T, const int& M, const float& tau_max,
     51                             const bool& reordering, std::vector<float>& tau) {
     52   // All cycles have the same fraction of the stopping time
     53   return fed_tau_by_cycle_time(T/(float)M,tau_max,reordering,tau);
     54 }
     55 
     56 //*************************************************************************************
     57 //*************************************************************************************
     58 
     59 /**
     60  * @brief This function allocates an array of the least number of time steps such
     61  * that a certain stopping time for the whole process can be obtained and fills it
     62  * it with the respective FED time step sizes for one cycle
     63  * The function returns the number of time steps per cycle or 0 on failure
     64  * @param t Desired cycle stopping time
     65  * @param tau_max Stability limit for the explicit scheme
     66  * @param reordering Reordering flag
     67  * @param tau The vector with the dynamic step sizes
     68  */
     69 int fed_tau_by_cycle_time(const float& t, const float& tau_max,
     70                           const bool& reordering, std::vector<float> &tau) {
     71   int n = 0;          // Number of time steps
     72   float scale = 0.0;  // Ratio of t we search to maximal t
     73 
     74   // Compute necessary number of time steps
     75   n = (int)(ceilf(sqrtf(3.0f*t/tau_max+0.25f)-0.5f-1.0e-8f)+ 0.5f);
     76   scale = 3.0f*t/(tau_max*(float)(n*(n+1)));
     77 
     78   // Call internal FED time step creation routine
     79   return fed_tau_internal(n,scale,tau_max,reordering,tau);
     80 }
     81 
     82 //*************************************************************************************
     83 //*************************************************************************************
     84 
     85 /**
     86  * @brief This function allocates an array of time steps and fills it with FED
     87  * time step sizes
     88  * The function returns the number of time steps per cycle or 0 on failure
     89  * @param n Number of internal steps
     90  * @param scale Ratio of t we search to maximal t
     91  * @param tau_max Stability limit for the explicit scheme
     92  * @param reordering Reordering flag
     93  * @param tau The vector with the dynamic step sizes
     94  */
     95 int fed_tau_internal(const int& n, const float& scale, const float& tau_max,
     96                      const bool& reordering, std::vector<float> &tau) {
     97   float c = 0.0, d = 0.0;     // Time savers
     98   vector<float> tauh;    // Helper vector for unsorted taus
     99 
    100   if (n <= 0) {
    101     return 0;
    102   }
    103 
    104   // Allocate memory for the time step size
    105   tau = vector<float>(n);
    106 
    107   if (reordering) {
    108     tauh = vector<float>(n);
    109   }
    110 
    111   // Compute time saver
    112   c = 1.0f / (4.0f * (float)n + 2.0f);
    113   d = scale * tau_max / 2.0f;
    114 
    115   // Set up originally ordered tau vector
    116   for (int k = 0; k < n; ++k) {
    117     float h = cosf((float)CV_PI * (2.0f * (float)k + 1.0f) * c);
    118 
    119     if (reordering) {
    120       tauh[k] = d / (h * h);
    121     }
    122     else {
    123       tau[k] = d / (h * h);
    124     }
    125   }
    126 
    127   // Permute list of time steps according to chosen reordering function
    128   int kappa = 0, prime = 0;
    129 
    130   if (reordering == true) {
    131     // Choose kappa cycle with k = n/2
    132     // This is a heuristic. We can use Leja ordering instead!!
    133     kappa = n / 2;
    134 
    135     // Get modulus for permutation
    136     prime = n + 1;
    137 
    138     while (!fed_is_prime_internal(prime)) {
    139       prime++;
    140     }
    141 
    142     // Perform permutation
    143     for (int k = 0, l = 0; l < n; ++k, ++l) {
    144       int index = 0;
    145       while ((index = ((k+1)*kappa) % prime - 1) >= n) {
    146         k++;
    147       }
    148 
    149       tau[l] = tauh[index];
    150     }
    151   }
    152 
    153   return n;
    154 }
    155 
    156 //*************************************************************************************
    157 //*************************************************************************************
    158 
    159 /**
    160  * @brief This function checks if a number is prime or not
    161  * @param number Number to check if it is prime or not
    162  * @return true if the number is prime
    163  */
    164 bool fed_is_prime_internal(const int& number) {
    165   bool is_prime = false;
    166 
    167   if (number <= 1) {
    168     return false;
    169   }
    170   else if (number == 1 || number == 2 || number == 3 || number == 5 || number == 7) {
    171     return true;
    172   }
    173   else if ((number % 2) == 0 || (number % 3) == 0 || (number % 5) == 0 || (number % 7) == 0) {
    174     return false;
    175   }
    176   else {
    177     is_prime = true;
    178     int upperLimit = (int)sqrt(1.0f + number);
    179     int divisor = 11;
    180 
    181     while (divisor <= upperLimit ) {
    182       if (number % divisor == 0)
    183       {
    184         is_prime = false;
    185       }
    186 
    187       divisor +=2;
    188     }
    189 
    190     return is_prime;
    191   }
    192 }
    193