1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
// BEGIN_COPYRIGHT
//
// Copyright (C) 1999 Allen Akin All Rights Reserved.
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the
// Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY
// KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
// WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
// PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL ALLEN AKIN BE
// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN
// AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
// OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// END_COPYRIGHT
// Image registration.
#include <cfloat>
#include "image.h"
#include <cmath> // for fabs
namespace GLEAN {
///////////////////////////////////////////////////////////////////////////////
// register: compare a reference image to the current (``test'') image.
//
// The reference image must be no larger than the current image, in
// both dimensions. Type doesn't matter, as both images will be
// converted to RGBA.
//
// The reference image will be slid into all possible positions over
// the current image, and the sum of the mean absolute errors for all
// four color channels computed at each position.
//
// Returns an Image::Registration struct that specifies the position at
// which the sum of mean absolute errors was minimal, plus the statistics
// at that position.
///////////////////////////////////////////////////////////////////////////////
Image::Registration
Image::reg(Image& ref) {
int wt = width(); // Width of test image, in pixels.
int ht = height(); // Height of test image, in pixels.
int wr = ref.width(); // Width of reference image, in pixels.
int hr = ref.height(); // Height of test image, in pixels.
int dh = ht - hr; // Difference in heights, in pixels.
int dw = wt - wr; // Difference in widths, in pixels.
int i;
if (dh < 0 || dw < 0)
throw RefImageTooLarge();
int wt4 = 4 * wt; // Width of test image, in RGBA samples.
int wr4 = 4 * wr; // Width of ref image, in RGBA samples.
int dw4 = 4 * dw; // Difference in widths, in samples.
double** testPix; // Buffers containing all the rows of
// the test image that need to be
// accessed concurrently.
// XXX sure would be nice to use auto_ptr to allocate this stuff,
// but it isn't supported in the STL that came with egcs 1.1.2.
// XXX testPix = new (double*) [dh + 1];
// VC 6 seems to misinterpret this as a c-style cast
testPix = new double* [dh + 1];
for (/*int */i = 0; i <= dh; ++i)
testPix[i] = new double [wt4];
double* refPix = new double [wr4];
// Buffer containing the one row of
// the reference image that's accessed
// at any given time.
BasicStats** stats; // Buffers containing a statistics-
// gathering structure for each of
// the possible reference image
// positions.
// XXX stats = new (BasicStats*) [dh + 1];
// VC 6 seems to misinterpret this as a c-style cast
stats = new BasicStats * [dh + 1];
for (/*int*/ i = 0; i <= dh; ++i)
stats[i] = new BasicStats [dw4 + 4];
// Prime the pump by unpacking the first few rows of the test image:
char* testRow = pixels();
for (/*int*/ i = 0; i < dh; ++i) {
unpack(wt, testPix[i], testRow);
testRow += rowSizeInBytes();
}
// Now accumulate statistics for one row of the reference image
// at a time, in all possible positions:
char* refRow = ref.pixels();
for (/*int*/ i = 0; i < hr; ++i) {
// Get the next row of the reference image:
ref.unpack(wr, refPix, refRow);
refRow += ref.rowSizeInBytes();
// Get the next row of the test image:
unpack(wt, testPix[dh], testRow);
testRow += rowSizeInBytes();
// Accumulate absolute error for R,G,B,A in all positions:
for (int j = 0; j <= dh; ++j)
for (int k = 0; k <= dw4; k += 4)
for (int m = 0; m < wr4; m += 4) {
stats[j][k+0].sample( fabs(
refPix[m+0]-testPix[j][m+k+0]));
stats[j][k+1].sample( fabs(
refPix[m+1]-testPix[j][m+k+1]));
stats[j][k+2].sample( fabs(
refPix[m+2]-testPix[j][m+k+2]));
stats[j][k+3].sample( fabs(
refPix[m+3]-testPix[j][m+k+3]));
}
}
// Now find the position for which the sum of the mean absolute errors
// is minimal:
double minErrorSum = DBL_MAX;
int minI = 0;
int minJ = 0;
for (/*int*/ i = 0; i <= dh; ++i)
for (int j = 0; j <= dw4; j += 4) {
double errorSum = stats[i][j+0].mean()
+ stats[i][j+1].mean()
+ stats[i][j+2].mean()
+ stats[i][j+3].mean();
if (errorSum < minErrorSum) {
minErrorSum = errorSum;
minI = i;
minJ = j;
}
}
Registration r;
r.wOffset = minJ / 4;
r.hOffset = minI;
r.stats[0] = stats[minI][minJ+0];
r.stats[1] = stats[minI][minJ+1];
r.stats[2] = stats[minI][minJ+2];
r.stats[3] = stats[minI][minJ+3];
// Clean up:
for (/*int*/ i = 0; i <= dh; ++i)
delete[] testPix[i];
delete[] testPix;
delete[] refPix;
for (/*int*/ i = 0; i <= dh; ++i)
delete[] stats[i];
delete[] stats;
return r;
} // Image::register
}; // namespace GLEAN
|
