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@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, Saarbrücken, 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  */
fed_tau_by_process_time(const float & T,const int & M,const float & tau_max,const bool & reordering,std::vector<float> & tau)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  */
fed_tau_by_cycle_time(const float & t,const float & tau_max,const bool & reordering,std::vector<float> & tau)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  */
fed_tau_internal(const int & n,const float & scale,const float & tau_max,const bool & reordering,std::vector<float> & tau)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  */
fed_is_prime_internal(const int & number)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