mnf.cpp
9.34 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
#include "linalg.h"
#include <iostream>
#include <stim/envi/envi.h>
#include "stim/envi/envi_header.h"
#include "stim/parser/arguments.h"
#include <stim/image/image.h>
#include <stim/math/matrix.h>
#include <time.h>
#include <fstream>
#include <thread>
#include "stim/envi/bil.h"
//LAPACKE support for Visual Studio
//#include <complex>
//#ifndef LAPACK_COMPLEX_CUSTOM
//#define LAPACK_COMPLEX_CUSTOM
//#define lapack_complex_float std::complex<float>
//#define lapack_complex_double std::complex<double>
//#endif
//#include "lapacke.h"
void progress_thread_envi(stim::envi* e); //progress bar threaded function
extern unsigned char* MASK; //pointer to an input file mask if one is provided
extern stim::envi ENVI; //input ENVI file to be processed
extern unsigned long long X, Y, B; //registers to quickly access the image size
///comparison function to compare two values corresponding to two memory addresses a and b
///is used in qsort function to sort elements of an array (ex.eigenvalues) and getting the corresponding indices.
int cmpfunc(const void * a, const void * b)
{
double *x = *(double **)a; double *y = *(double **)b;
return *x > *y ? 1 : -1;
}
/// multiply matrix A (nxm) by matrix B (mxp)
/// @param MatA is a pointer to the matrix A
/// @param MatB is a pointer to the matrix B
/// @param MatC is a pointer to pre-allocated memory of size [n * p * sizeof(double)] that stores the result of multiplication.
/// @param rowsX and @param columnsX are integers refering to the number of rows and columns of the matrix X respectively.
/// Note that columnA = rowsB
void multiply_2matrices(const double* MatA, int rowsA, int columnsA, const double* MatB, int rowsB, int columnsB, double* MatC, int rowsC, int columnsC) {
for (int i = 0; i < rowsA; i++) {
for (int j = 0; j < columnsB; j++) {
double sum = 0.0;
for (int k = 0; k < rowsB; k++)
sum = sum + MatA[i * columnsA + k] * MatB[k * columnsB + j];
MatC[i * columnsC + j] = sum; //resulting matrix C is [nxp]
}
}
}
void mnf(std::string outfile, int keptComponents, std::string NoiseFractions, int cuda_device){
///throw an Error in case of invalid values for the number of keptComponents
if (keptComponents <= 0 || keptComponents > B) {
std::cout << "ERROR: keptComponents must be in the range [0 " <<B<<"] (number of bands)"<< std::endl;
exit(1);
}
//throw an error if the ENVI file isn't BIP
if (ENVI.header.interleave != stim::envi_header::BIP) {
std::cout << "ERROR: MNF basis calculation requires a BIP file to make covariance calculation practical" << std::endl;
exit(1);
}
unsigned long long XY = X * Y; //number of pixels in the image
double* mu = (double*)malloc(sizeof(double) * B); //allocate space for the mean
std::thread t1(progress_thread_envi, &ENVI); //start the progress bar thread
std::cout << "Calculating mean spectrum..." << std::endl;
ENVI.mean_spectrum(mu, NULL, MASK, true); //calculate average spectrum
for (size_t b = 0; b < B; b++) { //for each band, test for non-finite values
if (!std::isfinite(mu[b])) //C++11 implementation
std::cout << "WARNING: the mean spectrum contains non-finite values in band " << b << ". Use --mask-finite to create a mask of finite values." << std::endl;
}
t1.join(); //wait for the progress bar thread to finish (it probably already is)
stim::matrix<double> cov(B, B); //allocate a double-precision covariance matrix
t1 = std::thread(progress_thread_envi, &ENVI); //start the progress bar thread
ENVI.co_matrix(cov.data(), mu, MASK, cuda_device, true); //calculate covariance matrix
t1.join(); //wait for the progress bar thread to finish (it probably already is)
stim::matrix<double> ncov(B, B); //allocate a double-precision noise covariance matrix
t1 = std::thread(progress_thread_envi, &ENVI); //start the progress bar thread
ENVI.coNoise_matrix(ncov.data(), mu, MASK, cuda_device, true); //calculate covariance of noise matrix
t1.join(); //wait for the progress bar thread to finish (it probably already is)
std::cout << std::endl << "Calculating the inverse covariance matrix S^(-1)...";
int *IPIV = (int*)malloc(sizeof(int) * B);
//LAPACKE_dgetrf(LAPACK_COL_MAJOR, (int)B, (int)B, cov.data(), (int)B, IPIV); //perform LU factorization
//LAPACKE_dgetri(LAPACK_COL_MAJOR, (int)B, cov.data(), (int)B, IPIV); //calculate matrix inverse
// REPLACED THE ABOVE LAPACKE functions with new CLAPACK functions in linalg.cpp
LINALG_dgetrf((int)B, (int)B, cov.data(), (int)B, IPIV); //perform LU factorization
LINALG_dgetri((int)B, cov.data(), (int)B, IPIV); //calculate matrix inverse
free(IPIV);
std::cout << "done." << std::endl;
std::cout << std::endl << "Calculating Q = S_n * S^(-1)..."; //calculate the Q matrix
stim::matrix<double> Q = ncov * cov;
std::cout << "done." << std::endl;
free(mu); //free memory for the mean spectrum
///computing the eigenvalues and left eigenvectors of matrix Qmat.
double *EigenvaluesReal = (double*)malloc(B * sizeof(double)); //allocate space for the real part of eigenvalues
double *EigenvaluesIm = (double*)malloc(B * sizeof(double)); //allocate space for the imaginary part of eigenvalues
std::cout << std::endl << "Calculating left eigenvectors X * A = v * X...";
stim::matrix<double> ev_left(B, B);
//LAPACKE_dgeev(LAPACK_COL_MAJOR, 'V', 'N', (int)B, Q.data(), (int)B, EigenvaluesReal, EigenvaluesIm, ev_left.data(), (int)B, 0, (int)B); //perform eigenvalue decomposition
// REPLACED THE ABOVE LAPACKE functions with new CLAPACK functions in linalg.cpp
LINALG_dgeev('V', 'N', (int)B, Q.data(), (int)B, EigenvaluesReal, EigenvaluesIm, ev_left.data(), (int)B, 0, (int)B); //perform eigenvalue decomposition
std::cout << "done." << std::endl;
std::cout << std::endl << "Sorting eigenvectors...";
///SORTING EIGENVECTORS CORRESPONDING TO THEIR EIGENVALUES
//EigenvaluesReal is the original array and "pointers" is the array of corresponding pointers, whose size is the size of the original array.
//Let cmp be the comparison function for qsort().
double **pointers = (double**)malloc(B * sizeof(double*));
for (int i = 0; i < B; i++) {
pointers[i] = &EigenvaluesReal[i]; // assign the address of integer.
}
///SORTING EIGENVALUES AND GETTING THE INDICES
qsort(pointers, B, sizeof(pointers[0]), cmpfunc);
size_t *EigenSortedIndices = (size_t*)malloc(B * sizeof(size_t)); //allocate space for the Indices of sorted eigenvalues
for (size_t i = 0; i < B; ++i)
EigenSortedIndices[i] = pointers[i] - EigenvaluesReal; //Calculate indices of the sorted eigenvalues
stim::matrix<double> As = ev_left.sort_cols(EigenSortedIndices);
stim::matrix<double> Ast = As.transpose();
std::cout << "done." << std::endl;
free(EigenvaluesReal);
free(EigenvaluesIm);
free(EigenSortedIndices);
std::cout << std::endl << "Removing noise components...";
///Noise removal and then backward transformation (builing one general/final basis matrix to be multiplied wih the raw data)
//General Transformation/Final Basis Matrix = gtm = (inv(As).') * (Keye) * (Ast) ;
//Note that the forward transform is done using transpose of mnf basis (Ast).However, for inverse transform, we use transpose of mnf basis inversed.
//creating a BxB matrix with the number of 'keptComponents' 1 on the diagonal and 0 elsewhere.
stim::matrix<double> K = stim::matrix<double>::I(B); //create an identity matrix
for (size_t i = keptComponents; i < B; i++) K(i, i) = 0.0; //assign zeros to components that will not be used
stim::matrix<double> K_Ast = K * Ast;
//calculate inverse of As matrix
int *IPIV2 = (int*)malloc(sizeof(int) * B);
//LAPACKE_dgetrf(LAPACK_COL_MAJOR, (int)B, (int)B, As.data(), (int)B, IPIV2);
//LAPACKE_dgetri(LAPACK_COL_MAJOR, (int)B, As.data(), (int)B, IPIV2);
// REPLACED THE ABOVE LAPACKE functions with new CLAPACK functions in linalg.cpp
LINALG_dgetrf((int)B, (int)B, As.data(), (int)B, IPIV2);
LINALG_dgetri((int)B, As.data(), (int)B, IPIV2);
free(IPIV2);
//calculate transpose of As inversed matrix (transpose of transformation matrix inveresed ) .
stim::matrix<double> Asit = As.transpose();
//multiplying As (inversed transposed matrix) by KeyeAst : (inv(As)).' * (KeyeAst) .
stim::matrix<double> G = Asit * K_Ast;
stim::matrix<double> G_t = G.transpose();
std::cout << "done." << std::endl;
std::ofstream basisfile(outfile);
basisfile << 0;
for (size_t i = 1; i < B; i++) basisfile << ", " << 0;
basisfile << std::endl;
basisfile << G.csv();
basisfile.close();
//output noise fractions
if (NoiseFractions.size() > 0) { //if a file name is specified for the noise fractions, save them
std::cout << std::endl << "Outputing noise fractions in a CSV file (infile-NoiseFractions)..." << std::endl;
//outputing the sorted Eigenvalues in a CSV file (Noise Fractions)
std::ofstream csv(NoiseFractions.c_str()); //open a CSV file to write the sorted Eigenvalues
csv << EigenSortedIndices[0] << ", " << EigenvaluesReal[EigenSortedIndices[0]]; //output the first variable
for (unsigned long long b = 1; b < B; b++)
csv << "\n" << EigenSortedIndices[b] << ", " << EigenvaluesReal[EigenSortedIndices[b]]; //output the next variable
csv.close();
}
} //end MNF