Subversion Repositories seema-scanner

Rev

Rev 236 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
184 jakw 1
#ifndef ALGORITHMTOOLS_H
2
#define ALGORITHMTOOLS_H
3
 
4
#include <opencv2/opencv.hpp>
236 jakw 5
#include <opencv2/core/ocl.hpp>
6
#include <QTime>
184 jakw 7
 
8
// Convert an unsigned binary number to reflected binary Gray code.
9
// Source: http://en.wikipedia.org/wiki/Gray_code
10
inline unsigned int binaryToGray(unsigned int num) {
11
    return (num >> 1) ^ num;
12
}
13
 
14
// Convert a reflected binary Gray code number to a binary number.
15
// Source: http://en.wikipedia.org/wiki/Gray_code
16
inline unsigned int grayToBinary(unsigned int num){
17
    unsigned int mask;
18
    for(mask = num >> 1; mask != 0; mask = mask >> 1)
19
        num = num ^ mask;
20
    return num;
21
}
22
 
23
// Return the Nth bit of an unsigned integer number
24
inline bool getBit(int decimal, int N){
25
 
26
    return decimal & 1 << (N-1);
27
}
28
 
29
// Return the position of the least significant bit that is set
30
inline int leastSignificantBitSet(int x){
31
  if(x == 0)
32
      return 0;
33
 
34
  int val = 1;
35
  while(x>>=1)
36
      val++;
37
 
38
  return val;
39
}
40
 
41
// Compute the power of a number (where the exponent is an integer)
42
inline unsigned int powi(int num, unsigned int exponent){
43
 
44
    if(exponent == 0)
45
        return 1;
46
 
47
    float res = num;
48
    for(unsigned int i=0; i<exponent-1; i++)
49
        res *= num;
50
 
51
    return res;
52
}
53
 
54
// Compute the power of a 2 (where the exponent is an integer)
55
inline unsigned int twopowi(unsigned int exponent){
56
    return 1 << exponent;
57
}
58
 
59
inline cv::Vec3b getColorSubpix(const cv::Mat& img, cv::Point2f pt){
60
    assert(!img.empty());
61
    assert(img.channels() == 3);
62
 
63
    int x = (int)pt.x;
64
    int y = (int)pt.y;
65
 
66
    int x0 = cv::borderInterpolate(x,   img.cols, cv::BORDER_REFLECT_101);
67
    int x1 = cv::borderInterpolate(x+1, img.cols, cv::BORDER_REFLECT_101);
68
    int y0 = cv::borderInterpolate(y,   img.rows, cv::BORDER_REFLECT_101);
69
    int y1 = cv::borderInterpolate(y+1, img.rows, cv::BORDER_REFLECT_101);
70
 
71
    float a = pt.x - (float)x;
72
    float c = pt.y - (float)y;
73
 
74
    uchar b = (uchar)cvRound((img.at<cv::Vec3b>(y0, x0)[0] * (1.f - a) + img.at<cv::Vec3b>(y0, x1)[0] * a) * (1.f - c)
75
                           + (img.at<cv::Vec3b>(y1, x0)[0] * (1.f - a) + img.at<cv::Vec3b>(y1, x1)[0] * a) * c);
76
    uchar g = (uchar)cvRound((img.at<cv::Vec3b>(y0, x0)[1] * (1.f - a) + img.at<cv::Vec3b>(y0, x1)[1] * a) * (1.f - c)
77
                           + (img.at<cv::Vec3b>(y1, x0)[1] * (1.f - a) + img.at<cv::Vec3b>(y1, x1)[1] * a) * c);
78
    uchar r = (uchar)cvRound((img.at<cv::Vec3b>(y0, x0)[2] * (1.f - a) + img.at<cv::Vec3b>(y0, x1)[2] * a) * (1.f - c)
79
                           + (img.at<cv::Vec3b>(y1, x0)[2] * (1.f - a) + img.at<cv::Vec3b>(y1, x1)[2] * a) * c);
80
 
81
    return cv::Vec3b(b, g, r);
82
}
83
 
84
inline cv::Mat computePhaseVector(unsigned int length, float phase, float pitch){
85
 
86
    cv::Mat phaseVector(length, 1, CV_8UC3);
87
    //phaseVector.setTo(0);
88
 
89
    const float pi = M_PI;
90
 
91
    // Loop through vector
92
    for(int i=0; i<phaseVector.rows; i++){
93
        // Amplitude of channels
94
        float amp = 0.5*(1+cos(2*pi*i/pitch - phase));
95
        phaseVector.at<cv::Vec3b>(i, 0) = cv::Vec3b(255.0*amp,255.0*amp,255.0*amp);
96
    }
97
 
98
    return phaseVector;
99
}
100
 
101
inline cv::Mat getPhase(const cv::Mat I1, const cv::Mat I2, const cv::Mat I3){
102
 
103
    cv::Mat_<float> I1_(I1);
104
    cv::Mat_<float> I2_(I2);
105
    cv::Mat_<float> I3_(I3);
106
 
107
    cv::Mat phase;
108
 
109
    // One call approach
110
    cv::phase(2.0*I1_-I3_-I2_, sqrt(3.0)*(I2_-I3_), phase);
111
    return phase;
112
 
113
}
114
 
115
// Phase unwrapping by means of a phase cue
207 flgw 116
inline cv::Mat unwrapWithCue(const cv::Mat &up, const cv::Mat &upCue, float nPhases){
184 jakw 117
 
118
    const float pi = M_PI;
119
 
187 jakw 120
    // Determine fringe order
184 jakw 121
    cv::Mat P = (upCue*nPhases-up)/(2.0*pi);
122
 
123
    // Round to integers
124
    P.convertTo(P, CV_8U);
125
    P.convertTo(P, CV_32F);
126
 
127
    // Add to phase
128
    cv::Mat upUnwrapped = up + P*2*pi;
129
 
130
    // Scale to range [0; 2pi]
131
    upUnwrapped *= 1.0/nPhases;
132
 
133
    return upUnwrapped;
134
}
135
 
136
// Absolute phase and magnitude from N frames
207 flgw 137
inline std::vector<cv::Mat> getDFTComponents(const std::vector<cv::Mat> &frames){
184 jakw 138
 
139
    unsigned int N = frames.size();
140
 
141
//    std::vector<cv::Mat> framesReverse = frames;
142
//    std::reverse(framesReverse.begin(), framesReverse.end());
143
 
144
    // DFT approach
145
    cv::Mat I;
146
    cv::merge(frames, I);
147
    unsigned int w = I.cols;
148
    unsigned int h = I.rows;
149
    I = I.reshape(1, h*w);
150
    I.convertTo(I, CV_32F);
236 jakw 151
 
245 jakw 152
    //QTime time;
153
    //time.start();
236 jakw 154
 
184 jakw 155
    cv::Mat fI;
156
    cv::dft(I, fI, cv::DFT_ROWS + cv::DFT_COMPLEX_OUTPUT);
236 jakw 157
 
245 jakw 158
    //std::cout << "elapsed: " << time.restart() << std::endl;
236 jakw 159
 
184 jakw 160
    fI = fI.reshape(N*2, h);
161
 
162
    std::vector<cv::Mat> fIcomp;
163
    cv::split(fI, fIcomp);
164
 
165
    return fIcomp;
166
 
167
}
168
 
245 jakw 169
// Determine mask based on positive monotony of phase in each row
170
inline void phaseMonotonyMask(cv::Mat& up, cv::Mat& mask){
171
 
172
    int nRows = up.rows;
173
    int nCols = up.cols;
174
 
175
    mask.create(nRows, nCols, CV_8UC1);
176
    mask.setTo(0);
177
 
178
    for(int r=0; r<nRows; r++){
179
 
180
        cv::Mat upRow = up.row(r);
181
        std::vector<int> b;
182
        //std::vector<float> lit;
183
        cv::Mat litMask(1, upRow.cols, CV_8UC1);
184
        litMask.setTo(0);
185
 
186
        // step 1: find longest increasing subsequence, (c) http://www.algorithmist.com
187
        std::vector<int> p(upRow.cols);
188
        int u, v;
189
 
190
        b.push_back(0);
191
 
192
        for(int i = 1; i < upRow.cols; i++){
193
 
194
            // If next element a[i] is greater than last element of
195
            // current longest subsequence, just push it at back of "b" and continue
196
            if (upRow.at<float>(0, b.back()) < upRow.at<float>(0, i)){
197
                p[i] = b.back();
198
                b.push_back(i);
199
                continue;
200
            }
201
 
202
            // Binary search to find the smallest element referenced by b which is just bigger than upRow[i]
203
            for(u = 0, v = b.size()-1; u < v;){
204
                int c = (u + v) / 2;
205
                if(upRow.at<float>(0, b[c]) < upRow.at<float>(0, i))
206
                    u = c+1;
207
                else
208
                    v = c;
209
            }
210
 
211
            // Update b if new value is smaller then previously referenced value
212
            if (upRow.at<float>(0, i) < upRow.at<float>(0, b[u])){
213
                if (u > 0) p[i] = b[u-1];
214
                b[u] = i;
215
            }
216
        }
217
 
218
        //lit.resize(b.size());
219
        for (u = b.size(), v = b.back(); u--; v = p[v]){
220
            litMask.at<uchar>(0, v) = 255;
221
            //lit[u] = v;
222
        }
223
 
224
        litMask.copyTo(mask.row(r));
225
    }
226
 
227
 
228
    return;
229
 
230
}
231
 
184 jakw 232
#endif // ALGORITHMTOOLS_H