cudapca.cpp
11.6 KB
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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* Copyright 2012 The MITRE Corporation *
* *
* Licensed under the Apache License, Version 2.0 (the "License"); *
* you may not use this file except in compliance with the License. *
* You may obtain a copy of the License at *
* *
* http://www.apache.org/licenses/LICENSE-2.0 *
* *
* Unless required by applicable law or agreed to in writing, software *
* distributed under the License is distributed on an "AS IS" BASIS, *
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
* See the License for the specific language governing permissions and *
* limitations under the License. *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#include <iostream>
using namespace std;
#include <unistd.h>
#include <QList>
#include <Eigen/Dense>
#include <opencv2/opencv.hpp>
using namespace cv;
#include <openbr/plugins/openbr_internal.h>
#include <openbr/core/common.h>
#include <openbr/core/eigenutils.h>
#include <openbr/core/opencvutils.h>
namespace br { namespace cuda { namespace pca {
void loadwrapper(float* evPtr, int evRows, int evCols, float* meanPtr, int meanElems);
void wrapper(void* src, void** dst);
}}}
namespace br
{
/*!
* \ingroup transforms
* \brief Projects input into learned Principal Component Analysis subspace using CUDA.
* \author Brendan Klare \cite bklare
* \author Josh Klontz \cite jklontz
* \author Colin Heinzmann \cite DepthDeluxe
*
* \br_property float keep Options are: [keep < 0 - All eigenvalues are retained, keep == 0 - No PCA is performed and the eigenvectors form an identity matrix, 0 < keep < 1 - Keep is the fraction of the variance to retain, keep >= 1 - keep is the number of leading eigenvectors to retain] Default is 0.95.
* \br_property int drop The number of leading eigen-dimensions to drop.
* \br_property bool whiten Whether or not to perform PCA whitening (i.e., normalize variance of each dimension to unit norm)
*/
class CUDAPCATransform : public Transform
{
Q_OBJECT
protected:
Q_PROPERTY(float keep READ get_keep WRITE set_keep RESET reset_keep STORED false)
Q_PROPERTY(int drop READ get_drop WRITE set_drop RESET reset_drop STORED false)
Q_PROPERTY(bool whiten READ get_whiten WRITE set_whiten RESET reset_whiten STORED false)
BR_PROPERTY(float, keep, 0.95)
BR_PROPERTY(int, drop, 0)
BR_PROPERTY(bool, whiten, false)
Eigen::VectorXf mean, eVals;
Eigen::MatrixXf eVecs;
int originalRows;
public:
CUDAPCATransform() : keep(0.95), drop(0), whiten(false) {}
private:
double residualReconstructionError(const Template &src) const
{
Template proj;
project(src, proj);
Eigen::Map<const Eigen::VectorXf> srcMap(src.m().ptr<float>(), src.m().rows*src.m().cols);
Eigen::Map<Eigen::VectorXf> projMap(proj.m().ptr<float>(), keep);
return (srcMap - mean).squaredNorm() - projMap.squaredNorm();
}
void train(const TemplateList &cudaTrainingSet)
{
// copy the data back from the graphics card so the training can be done on the CPU
const int instances = cudaTrainingSet.size(); // get the number of training set instances
QList<Template> trainingQlist;
for(int i=0; i<instances; i++) {
Template currentTemplate = cudaTrainingSet[i];
void* const* srcDataPtr = currentTemplate.m().ptr<void*>();
void* cudaMemPtr = srcDataPtr[0];
int rows = *((int*)srcDataPtr[1]);
int cols = *((int*)srcDataPtr[2]);
int type = *((int*)srcDataPtr[3]);
Mat mat = Mat(rows, cols, type);
trainingQlist.append(Template(mat));
}
// assemble a TemplateList from the list of data
TemplateList trainingSet(trainingQlist);
if (trainingSet.first().m().type() != CV_32FC1) {
qFatal("Requires single channel 32-bit floating point matrices.");
}
originalRows = trainingSet.first().m().rows; // get number of rows of first image
int dimsIn = trainingSet.first().m().rows * trainingSet.first().m().cols; // get the size of the first image
// Map into 64-bit Eigen matrix
Eigen::MatrixXd data(dimsIn, instances); // create a mat
for (int i=0; i<instances; i++) {
data.col(i) = Eigen::Map<const Eigen::MatrixXf>(trainingSet[i].m().ptr<float>(), dimsIn, 1).cast<double>();
}
trainCore(data);
}
void project(const Template &src, Template &dst) const
{
void* const* srcDataPtr = src.m().ptr<void*>();
int rows = *((int*)srcDataPtr[1]);
int cols = *((int*)srcDataPtr[2]);
int type = *((int*)srcDataPtr[3]);
if (type != CV_32FC1) {
cout << "ERR: Invalid image type" << endl;
return;
}
Mat dstMat = Mat(src.m().rows, src.m().cols, src.m().type());
void** dstDataPtr = dstMat.ptr<void*>();
dstDataPtr[1] = srcDataPtr[1]; *((int*)dstDataPtr[1]) = 1;
dstDataPtr[2] = srcDataPtr[2]; *((int*)dstDataPtr[2]) = keep;
dstDataPtr[3] = srcDataPtr[3];
br::cuda::pca::wrapper(srcDataPtr[0], &dstDataPtr[0]);
dst = dstMat;
//dst = cv::Mat(1, keep, CV_32FC1);
// perform the operation on the graphics card
//cuda::cudapca_projectwrapper((float*)src.m().ptr<float>(), (float*)dst.m().ptr<float>());
// Map Eigen into OpenCV
//Mat cpuDst = cv::Mat(1, keep, CV_32FC1);
//Eigen::Map<const Eigen::MatrixXf> inMap(src.m().ptr<float>(), src.m().rows*src.m().cols, 1);
//Eigen::Map<Eigen::MatrixXf> outMap(cpuDst.ptr<float>(), keep, 1);
// Do projection
//outMap = eVecs.transpose() * (inMap - mean);
}
void store(QDataStream &stream) const
{
stream << keep << drop << whiten << originalRows << mean << eVals << eVecs;
}
void load(QDataStream &stream)
{
Eigen::MatrixXf originalEVecs;
stream >> keep >> drop >> whiten >> originalRows >> mean >> eVals >> originalEVecs;
// perform transpose before copying over
eVecs = originalEVecs; //originalEVecs.transpose();
cout << "Mean Dimensions" << endl;
cout << "\tRows: " << mean.rows() << " Cols: " << mean.cols() << endl;
cout << "eVecs Dimensions" << endl;
cout << "\tRows: " << eVecs.rows() << " Cols: " << eVecs.cols() << endl;
cout << "eVals Dimensions" << endl;
cout << "\tRows: " << eVals.rows() << " Cols: " << eVals.cols() << endl;
cout << "Keep: " << keep << endl;
cout << "Mean first value: " << mean(0, 0) << endl;
// TODO(colin): use Eigen Map class to generate map files so we don't have to copy the data
// serialize the eigenvectors
float* evBuffer = new float[eVecs.rows() * eVecs.cols()];
for (int i=0; i < eVecs.rows(); i++) {
for (int j=0; j < eVecs.cols(); j++) {
evBuffer[i*eVecs.cols() + j] = eVecs(i, j);
}
}
// serialize the mean
float* meanBuffer = new float[mean.rows() * mean.cols()];
for (int i=0; i < mean.rows(); i++) {
for (int j=0; j < mean.cols(); j++) {
meanBuffer[i*mean.cols() + j] = mean(i, j);
}
}
// call the wrapper function
br::cuda::pca::loadwrapper(evBuffer, eVecs.rows(), eVecs.cols(), meanBuffer, mean.rows()*mean.cols());
delete evBuffer;
delete meanBuffer;
}
protected:
void trainCore(Eigen::MatrixXd data)
{
int dimsIn = data.rows();
int instances = data.cols();
const bool dominantEigenEstimation = (dimsIn > instances);
Eigen::MatrixXd allEVals, allEVecs;
if (keep != 0) {
// Compute and remove mean
mean = Eigen::VectorXf(dimsIn);
for (int i=0; i<dimsIn; i++) mean(i) = data.row(i).sum() / (float)instances;
for (int i=0; i<dimsIn; i++) data.row(i).array() -= mean(i);
// Calculate covariance matrix
Eigen::MatrixXd cov;
if (dominantEigenEstimation) cov = data.transpose() * data / (instances-1.0);
else cov = data * data.transpose() / (instances-1.0);
// Compute eigendecomposition. Returns eigenvectors/eigenvalues in increasing order by eigenvalue.
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eSolver(cov);
allEVals = eSolver.eigenvalues();
allEVecs = eSolver.eigenvectors();
if (dominantEigenEstimation) allEVecs = data * allEVecs;
} else {
// Null case
mean = Eigen::VectorXf::Zero(dimsIn);
allEVecs = Eigen::MatrixXd::Identity(dimsIn, dimsIn);
allEVals = Eigen::VectorXd::Ones(dimsIn);
}
if (keep <= 0) {
keep = dimsIn - drop;
} else if (keep < 1) {
// Keep eigenvectors that retain a certain energy percentage.
const double totalEnergy = allEVals.sum();
if (totalEnergy == 0) {
keep = 0;
} else {
double currentEnergy = 0;
int i=0;
while ((currentEnergy / totalEnergy < keep) && (i < allEVals.rows())) {
currentEnergy += allEVals(allEVals.rows()-(i+1));
i++;
}
keep = i - drop;
}
} else {
if (keep + drop > allEVals.rows()) {
qWarning("Insufficient samples, needed at least %d but only got %d.", (int)keep + drop, (int)allEVals.rows());
keep = allEVals.rows() - drop;
}
}
// Keep highest energy vectors
eVals = Eigen::VectorXf((int)keep, 1);
eVecs = Eigen::MatrixXf(allEVecs.rows(), (int)keep);
for (int i=0; i<keep; i++) {
int index = allEVals.rows()-(i+drop+1);
eVals(i) = allEVals(index);
eVecs.col(i) = allEVecs.col(index).cast<float>() / allEVecs.col(index).norm();
if (whiten) eVecs.col(i) /= sqrt(eVals(i));
}
// Debug output
if (Globals->verbose) qDebug() << "PCA Training:\n\tDimsIn =" << dimsIn << "\n\tKeep =" << keep;
}
void writeEigenVectors(const Eigen::MatrixXd &allEVals, const Eigen::MatrixXd &allEVecs) const
{
const int originalCols = mean.rows() / originalRows;
{ // Write out mean image
cv::Mat out(originalRows, originalCols, CV_32FC1);
Eigen::Map<Eigen::MatrixXf> outMap(out.ptr<float>(), mean.rows(), 1);
outMap = mean.col(0);
// OpenCVUtils::saveImage(out, Globals->Debug+"/PCA/eigenVectors/mean.png");
}
// Write out sample eigen vectors (16 highest, 8 lowest), filename = eigenvalue.
for (int k=0; k<(int)allEVals.size(); k++) {
if ((k < 8) || (k >= (int)allEVals.size()-16)) {
cv::Mat out(originalRows, originalCols, CV_64FC1);
Eigen::Map<Eigen::MatrixXd> outMap(out.ptr<double>(), mean.rows(), 1);
outMap = allEVecs.col(k);
// OpenCVUtils::saveImage(out, Globals->Debug+"/PCA/eigenVectors/"+QString::number(allEVals(k),'f',0)+".png");
}
}
}
};
BR_REGISTER(Transform, CUDAPCATransform)
} // namespace br
#include "cuda/cudapca.moc"