1/*
2* This sample demonstrates the use of the function
3* findTransformECC that implements the image alignment ECC algorithm
4*
5*
6* The demo loads an image (defaults to ../data/fruits.jpg) and it artificially creates
7* a template image based on the given motion type. When two images are given,
8* the first image is the input image and the second one defines the template image.
9* In the latter case, you can also parse the warp's initialization.
10*
11* Input and output warp files consist of the raw warp (transform) elements
12*
13* Authors: G. Evangelidis, INRIA, Grenoble, France
14*          M. Asbach, Fraunhofer IAIS, St. Augustin, Germany
15*/
16#include <opencv2/imgcodecs.hpp>
17#include <opencv2/highgui.hpp>
18#include <opencv2/video.hpp>
19#include <opencv2/imgproc.hpp>
20#include <opencv2/core/utility.hpp>
21
22#include <stdio.h>
23#include <string>
24#include <time.h>
25#include <iostream>
26#include <fstream>
27
28using namespace cv;
29using namespace std;
30
31static void help(void);
32static int readWarp(string iFilename, Mat& warp, int motionType);
33static int saveWarp(string fileName, const Mat& warp, int motionType);
34static void draw_warped_roi(Mat& image, const int width, const int height, Mat& W);
35
36#define HOMO_VECTOR(H, x, y)\
37    H.at<float>(0,0) = (float)(x);\
38    H.at<float>(1,0) = (float)(y);\
39    H.at<float>(2,0) = 1.;
40
41#define GET_HOMO_VALUES(X, x, y)\
42    (x) = static_cast<float> (X.at<float>(0,0)/X.at<float>(2,0));\
43    (y) = static_cast<float> (X.at<float>(1,0)/X.at<float>(2,0));
44
45
46const std::string keys =
47    "{@inputImage    | ../data/fruits.jpg | input image filename }"
48    "{@templateImage |               | template image filename (optional)}"
49    "{@inputWarp     |               | input warp (matrix) filename (optional)}"
50    "{n numOfIter    | 50            | ECC's iterations }"
51    "{e epsilon      | 0.0001        | ECC's convergence epsilon }"
52    "{o outputWarp   | outWarp.ecc   | output warp (matrix) filename }"
53    "{m motionType   | affine        | type of motion (translation, euclidean, affine, homography) }"
54    "{v verbose      | 0             | display initial and final images }"
55    "{w warpedImfile | warpedECC.png | warped input image }"
56;
57
58
59static void help(void)
60{
61
62    cout << "\nThis file demostrates the use of the ECC image alignment algorithm. When one image"
63        " is given, the template image is artificially formed by a random warp. When both images"
64        " are given, the initialization of the warp by command line parsing is possible. "
65        "If inputWarp is missing, the identity transformation initializes the algorithm. \n" << endl;
66
67    cout << "\nUsage example (one image): \n./ecc ../data/fruits.jpg -o=outWarp.ecc "
68        "-m=euclidean -e=1e-6 -N=70 -v=1 \n" << endl;
69
70    cout << "\nUsage example (two images with initialization): \n./ecc yourInput.png yourTemplate.png "
71        "yourInitialWarp.ecc -o=outWarp.ecc -m=homography -e=1e-6 -N=70 -v=1 -w=yourFinalImage.png \n" << endl;
72
73}
74
75static int readWarp(string iFilename, Mat& warp, int motionType){
76
77    // it reads from file a specific number of raw values:
78    // 9 values for homography, 6 otherwise
79    CV_Assert(warp.type()==CV_32FC1);
80    int numOfElements;
81    if (motionType==MOTION_HOMOGRAPHY)
82        numOfElements=9;
83    else
84        numOfElements=6;
85
86    int i;
87    int ret_value;
88
89    ifstream myfile(iFilename.c_str());
90    if (myfile.is_open()){
91        float* matPtr = warp.ptr<float>(0);
92        for(i=0; i<numOfElements; i++){
93            myfile >> matPtr[i];
94        }
95        ret_value = 1;
96    }
97    else {
98        cout << "Unable to open file " << iFilename.c_str() << endl;
99        ret_value = 0;
100    }
101    return ret_value;
102}
103
104static int saveWarp(string fileName, const Mat& warp, int motionType)
105{
106    // it saves the raw matrix elements in a file
107    CV_Assert(warp.type()==CV_32FC1);
108
109    const float* matPtr = warp.ptr<float>(0);
110    int ret_value;
111
112    ofstream outfile(fileName.c_str());
113    if( !outfile ) {
114        cerr << "error in saving "
115            << "Couldn't open file '" << fileName.c_str() << "'!" << endl;
116        ret_value = 0;
117    }
118    else {//save the warp's elements
119        outfile << matPtr[0] << " " << matPtr[1] << " " << matPtr[2] << endl;
120        outfile << matPtr[3] << " " << matPtr[4] << " " << matPtr[5] << endl;
121        if (motionType==MOTION_HOMOGRAPHY){
122            outfile << matPtr[6] << " " << matPtr[7] << " " << matPtr[8] << endl;
123        }
124        ret_value = 1;
125    }
126    return ret_value;
127
128}
129
130
131static void draw_warped_roi(Mat& image, const int width, const int height, Mat& W)
132{
133    Point2f top_left, top_right, bottom_left, bottom_right;
134
135    Mat  H = Mat (3, 1, CV_32F);
136    Mat  U = Mat (3, 1, CV_32F);
137
138    Mat warp_mat = Mat::eye (3, 3, CV_32F);
139
140    for (int y = 0; y < W.rows; y++)
141        for (int x = 0; x < W.cols; x++)
142            warp_mat.at<float>(y,x) = W.at<float>(y,x);
143
144    //warp the corners of rectangle
145
146    // top-left
147    HOMO_VECTOR(H, 1, 1);
148    gemm(warp_mat, H, 1, 0, 0, U);
149    GET_HOMO_VALUES(U, top_left.x, top_left.y);
150
151    // top-right
152    HOMO_VECTOR(H, width, 1);
153    gemm(warp_mat, H, 1, 0, 0, U);
154    GET_HOMO_VALUES(U, top_right.x, top_right.y);
155
156    // bottom-left
157    HOMO_VECTOR(H, 1, height);
158    gemm(warp_mat, H, 1, 0, 0, U);
159    GET_HOMO_VALUES(U, bottom_left.x, bottom_left.y);
160
161    // bottom-right
162    HOMO_VECTOR(H, width, height);
163    gemm(warp_mat, H, 1, 0, 0, U);
164    GET_HOMO_VALUES(U, bottom_right.x, bottom_right.y);
165
166    // draw the warped perimeter
167    line(image, top_left, top_right, Scalar(255,0,255));
168    line(image, top_right, bottom_right, Scalar(255,0,255));
169    line(image, bottom_right, bottom_left, Scalar(255,0,255));
170    line(image, bottom_left, top_left, Scalar(255,0,255));
171}
172
173int main (const int argc, const char * argv[])
174{
175
176    CommandLineParser parser(argc, argv, keys);
177    parser.about("ECC demo");
178
179    if (argc<2) {
180        parser.printMessage();
181        help();
182        return 1;
183    }
184
185    string imgFile = parser.get<string>(0);
186    string tempImgFile = parser.get<string>(1);
187    string inWarpFile = parser.get<string>(2);
188
189    int number_of_iterations = parser.get<int>("n");
190    double termination_eps = parser.get<double>("e");
191    string warpType = parser.get<string>("m");
192    int verbose = parser.get<int>("v");
193    string finalWarp = parser.get<string>("o");
194    string warpedImFile = parser.get<string>("w");
195
196    if (!(warpType == "translation" || warpType == "euclidean"
197        || warpType == "affine" || warpType == "homography"))
198    {
199        cerr << "Invalid motion transformation" << endl;
200        return -1;
201    }
202
203    int mode_temp;
204    if (warpType == "translation")
205        mode_temp = MOTION_TRANSLATION;
206    else if (warpType == "euclidean")
207        mode_temp = MOTION_EUCLIDEAN;
208    else if (warpType == "affine")
209        mode_temp = MOTION_AFFINE;
210    else
211        mode_temp = MOTION_HOMOGRAPHY;
212
213    Mat inputImage = imread(imgFile,0);
214    if (inputImage.empty())
215    {
216        cerr << "Unable to load the inputImage" <<  endl;
217        return -1;
218    }
219
220    Mat target_image;
221    Mat template_image;
222
223    if (tempImgFile!="") {
224        inputImage.copyTo(target_image);
225        template_image = imread(tempImgFile,0);
226        if (template_image.empty()){
227            cerr << "Unable to load the template image" << endl;
228            return -1;
229        }
230
231    }
232    else{ //apply random waro to input image
233        resize(inputImage, target_image, Size(216, 216));
234        Mat warpGround;
235        cv::RNG rng;
236        double angle;
237        switch (mode_temp) {
238        case MOTION_TRANSLATION:
239            warpGround = (Mat_<float>(2,3) << 1, 0, (rng.uniform(10.f, 20.f)),
240                0, 1, (rng.uniform(10.f, 20.f)));
241            warpAffine(target_image, template_image, warpGround,
242                Size(200,200), INTER_LINEAR + WARP_INVERSE_MAP);
243            break;
244        case MOTION_EUCLIDEAN:
245            angle = CV_PI/30 + CV_PI*rng.uniform((double)-2.f, (double)2.f)/180;
246
247            warpGround = (Mat_<float>(2,3) << cos(angle), -sin(angle), (rng.uniform(10.f, 20.f)),
248                sin(angle), cos(angle), (rng.uniform(10.f, 20.f)));
249            warpAffine(target_image, template_image, warpGround,
250                Size(200,200), INTER_LINEAR + WARP_INVERSE_MAP);
251            break;
252        case MOTION_AFFINE:
253
254            warpGround = (Mat_<float>(2,3) << (1-rng.uniform(-0.05f, 0.05f)),
255                (rng.uniform(-0.03f, 0.03f)), (rng.uniform(10.f, 20.f)),
256                (rng.uniform(-0.03f, 0.03f)), (1-rng.uniform(-0.05f, 0.05f)),
257                (rng.uniform(10.f, 20.f)));
258            warpAffine(target_image, template_image, warpGround,
259                Size(200,200), INTER_LINEAR + WARP_INVERSE_MAP);
260            break;
261        case MOTION_HOMOGRAPHY:
262            warpGround = (Mat_<float>(3,3) << (1-rng.uniform(-0.05f, 0.05f)),
263                (rng.uniform(-0.03f, 0.03f)), (rng.uniform(10.f, 20.f)),
264                (rng.uniform(-0.03f, 0.03f)), (1-rng.uniform(-0.05f, 0.05f)),(rng.uniform(10.f, 20.f)),
265                (rng.uniform(0.0001f, 0.0003f)), (rng.uniform(0.0001f, 0.0003f)), 1.f);
266            warpPerspective(target_image, template_image, warpGround,
267                Size(200,200), INTER_LINEAR + WARP_INVERSE_MAP);
268            break;
269        }
270    }
271
272
273    const int warp_mode = mode_temp;
274
275    // initialize or load the warp matrix
276    Mat warp_matrix;
277    if (warpType == "homography")
278        warp_matrix = Mat::eye(3, 3, CV_32F);
279    else
280        warp_matrix = Mat::eye(2, 3, CV_32F);
281
282    if (inWarpFile!=""){
283        int readflag = readWarp(inWarpFile, warp_matrix, warp_mode);
284        if ((!readflag) || warp_matrix.empty())
285        {
286            cerr << "-> Check warp initialization file" << endl << flush;
287            return -1;
288        }
289    }
290    else {
291
292        printf("\n ->Perfomarnce Warning: Identity warp ideally assumes images of "
293            "similar size. If the deformation is strong, the identity warp may not "
294            "be a good initialization. \n");
295
296    }
297
298    if (number_of_iterations > 200)
299        cout << "-> Warning: too many iterations " << endl;
300
301    if (warp_mode != MOTION_HOMOGRAPHY)
302        warp_matrix.rows = 2;
303
304    // start timing
305    const double tic_init = (double) getTickCount ();
306    double cc = findTransformECC (template_image, target_image, warp_matrix, warp_mode,
307        TermCriteria (TermCriteria::COUNT+TermCriteria::EPS,
308        number_of_iterations, termination_eps));
309
310    if (cc == -1)
311    {
312        cerr << "The execution was interrupted. The correlation value is going to be minimized." << endl;
313        cerr << "Check the warp initialization and/or the size of images." << endl << flush;
314    }
315
316    // end timing
317    const double toc_final  = (double) getTickCount ();
318    const double total_time = (toc_final-tic_init)/(getTickFrequency());
319    if (verbose){
320        cout << "Alignment time (" << warpType << " transformation): "
321            << total_time << " sec" << endl << flush;
322        //  cout << "Final correlation: " << cc << endl << flush;
323
324    }
325
326    // save the final warp matrix
327    saveWarp(finalWarp, warp_matrix, warp_mode);
328
329    if (verbose){
330        cout << "\nThe final warp has been saved in the file: " << finalWarp << endl << flush;
331    }
332
333    // save the final warped image
334    Mat warped_image = Mat(template_image.rows, template_image.cols, CV_32FC1);
335    if (warp_mode != MOTION_HOMOGRAPHY)
336        warpAffine      (target_image, warped_image, warp_matrix, warped_image.size(),
337        INTER_LINEAR + WARP_INVERSE_MAP);
338    else
339        warpPerspective (target_image, warped_image, warp_matrix, warped_image.size(),
340        INTER_LINEAR + WARP_INVERSE_MAP);
341
342    //save the warped image
343    imwrite(warpedImFile, warped_image);
344
345    // display resulting images
346    if (verbose)
347    {
348
349        cout << "The warped image has been saved in the file: " << warpedImFile << endl << flush;
350
351        namedWindow ("image",    WINDOW_AUTOSIZE);
352        namedWindow ("template", WINDOW_AUTOSIZE);
353        namedWindow ("warped image",   WINDOW_AUTOSIZE);
354        namedWindow ("error (black: no error)", WINDOW_AUTOSIZE);
355
356        moveWindow  ("template", 350, 350);
357        moveWindow  ("warped image",   600, 300);
358        moveWindow  ("error (black: no error)", 900, 300);
359
360        // draw boundaries of corresponding regions
361        Mat identity_matrix = Mat::eye(3,3,CV_32F);
362
363        draw_warped_roi (target_image,   template_image.cols-2, template_image.rows-2, warp_matrix);
364        draw_warped_roi (template_image, template_image.cols-2, template_image.rows-2, identity_matrix);
365
366        Mat errorImage;
367        subtract(template_image, warped_image, errorImage);
368        double max_of_error;
369        minMaxLoc(errorImage, NULL, &max_of_error);
370
371        // show images
372        cout << "Press any key to exit the demo (you might need to click on the images before)." << endl << flush;
373
374        imshow ("image",    target_image);
375        waitKey (200);
376        imshow ("template", template_image);
377        waitKey (200);
378        imshow ("warped image",   warped_image);
379        waitKey(200);
380        imshow ("error (black: no error)",  abs(errorImage)*255/max_of_error);
381        waitKey(0);
382
383    }
384
385    // done
386    return 0;
387}
388