Skip to content

Commit d1dbf4b

Browse files
jcollfontjcollfont
committed
updates with Fernando. Fixed many bugs in the code. currently not compiling
1 parent 386ac6f commit d1dbf4b

File tree

9 files changed

+180
-170
lines changed

9 files changed

+180
-170
lines changed

‎src/Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovImpl_child.cc‎

Lines changed: 33 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,7 @@ using namespace SCIRun::Core::Algorithms::Inverse;
5757

5858
/////////////////////////
5959
///////// compute Inverse solution
60-
DenseColumnMatrix SolveInverseProblemWithTikhonovImpl_child::computeInverseSolution( double lambda_sq, bool inverseCalculation)
60+
DenseMatrix SolveInverseProblemWithTikhonovImpl_child::computeInverseSolution( double lambda_sq, bool inverseCalculation) const
6161
{
6262
//............................
6363
// OPERATIONS PERFORMED IN THIS SECTION:
@@ -73,10 +73,11 @@ using namespace SCIRun::Core::Algorithms::Inverse;
7373

7474
const int sizeB = M1.ncols();
7575
const int sizeSolution = M3.nrows();
76+
const int numTimeSamples = y.ncols();
7677
DenseMatrix inverseG(sizeB,sizeB);
7778

7879
DenseColumnMatrix b(sizeB);
79-
DenseColumnMatrix solution(sizeSolution);
80+
DenseMatrix solution(sizeSolution,numTimeSamples);
8081
DenseMatrix G;
8182

8283
G = M1 + lambda_sq * M2;
@@ -98,16 +99,23 @@ using namespace SCIRun::Core::Algorithms::Inverse;
9899

99100
/////// precomputeInverseMatrices
100101
///////////////
101-
void SolveInverseProblemWithTikhonovImpl_child::preAlocateInverseMatrices(SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_)
102+
void SolveInverseProblemWithTikhonovImpl_child::preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_)
102103
{
103104

104105
// TODO: use DimensionMismatch exception where appropriate
105106
// DIMENSION CHECK!!
106-
const int M = forwardMatrix_->nrows();
107-
const int N = forwardMatrix_->ncols();
107+
const int M = forwardMatrix_.nrows();
108+
const int N = forwardMatrix_.ncols();
109+
110+
int x = 0;
108111

109112
// PREALOCATE VARIABLES and MATRICES
110-
DenseMatrix forward_transpose = forwardMatrix_->transpose();
113+
DenseMatrix forward_transpose = forwardMatrix_.transpose();
114+
115+
// get Parameters
116+
auto regularizationChoice_ = get(regularizationChoice).toInt();
117+
auto regularizationSolutionSubcase_ = get(regularizationSolutionSubcase).toInt();
118+
auto regularizationResidualSubcase_ = get(regularizationResidualSubcase).toInt();
111119

112120
// select underdetermined case if user decides so or the option is set to automatic and number of measurements is smaller than number of unknowns.
113121
if ( ( (M < N) && (regularizationChoice_ == automatic) ) || (regularizationChoice_ == underdetermined))
@@ -132,7 +140,7 @@ using namespace SCIRun::Core::Algorithms::Inverse;
132140

133141
// DEFINITIONS AND PREALOCATION OF SOURCE REGULARIZATION MATRIX 'R'
134142
// if R does not exist, set as identity of size equal to N (columns of fwd matrix)
135-
if (!sourceWeighting_)
143+
if (true)//(&sourceWeighting_==NULL)
136144
{
137145
RRtr = DenseMatrix::Identity(N, N);
138146
iRRtr = RRtr;
@@ -143,28 +151,19 @@ using namespace SCIRun::Core::Algorithms::Inverse;
143151
// if provided the non-squared version of R
144152
if( regularizationSolutionSubcase_==solution_constrained )
145153
{
146-
RRtr = sourceWeighting_->transpose() * *sourceWeighting_;
154+
RRtr = sourceWeighting_.transpose() * sourceWeighting_;
147155
}
148156
// otherwise, if the source regularization is provided as the squared version (RR^T)
149157
else if ( regularizationSolutionSubcase_==solution_constrained_squared )
150158
{
151-
RRtr = *sourceWeighting_;
159+
RRtr = sourceWeighting_;
152160
}
153161

154162
// check if squared regularization matrix is invertible
155163
if ( !RRtr.fullPivLu().isInvertible() )
156164
{
157-
const std::string errorMessage("Regularization matrix in the source space is not invertible.");
158-
if (pr_)
159-
{
160-
pr_->error(errorMessage);
161-
}
162-
else
163-
{
164-
std::cerr << errorMessage << std::endl;
165-
}
166165

167-
THROW_ALGORITHM_INPUT_ERROR_SIMPLE(errorMessage);
166+
THROW_ALGORITHM_INPUT_ERROR_SIMPLE("Regularization matrix in the source space is not invertible.");
168167
}
169168

170169
// COMPUTE inverse
@@ -175,7 +174,7 @@ using namespace SCIRun::Core::Algorithms::Inverse;
175174

176175
// DEFINITIONS AND PREALOCATIONS OF MEASUREMENTS COVARIANCE MATRIX 'C'
177176
// if C does not exist, set as identity of size equal to M (rows of fwd matrix)
178-
if (!sensorWeighting_)
177+
if (true)//(&sensorWeighting_==NULL)
179178
{
180179
CCtr = DenseMatrix::Identity(M, M);
181180
iCCtr = CCtr;
@@ -187,30 +186,21 @@ using namespace SCIRun::Core::Algorithms::Inverse;
187186
if (regularizationResidualSubcase_ == residual_constrained)
188187
{
189188
// check that the matrix is of appropriate size (equal number of rows as rows in fwd matrix)
190-
if(M != sensorWeighting_->ncols())
189+
if(M != sensorWeighting_.ncols())
191190
{
192-
CCtr = sensorWeighting_->transpose() * *sensorWeighting_;
191+
CCtr = sensorWeighting_.transpose() * sensorWeighting_;
193192
}
194193
}
195194
// otherwise if the source covariance matrix is provided in squared form
196195
else if ( regularizationResidualSubcase_ == residual_constrained_squared )
197196
{
198-
CCtr = *sensorWeighting_;
197+
CCtr = sensorWeighting_;
199198
}
200199

201200
// check if squared regularization matrix is invertible
202201
if ( !CCtr.fullPivLu().isInvertible() )
203202
{
204-
const std::string errorMessage("Residual covariance matrix is not invertible.");
205-
if (pr_)
206-
{
207-
pr_->error(errorMessage);
208-
}
209-
else
210-
{
211-
std::cerr << errorMessage << std::endl;
212-
}
213-
THROW_ALGORITHM_INPUT_ERROR_SIMPLE(errorMessage);
203+
THROW_ALGORITHM_INPUT_ERROR_SIMPLE("Residual covariance matrix is not invertible.");
214204
}
215205
iCCtr = CCtr.inverse().eval();
216206

@@ -220,7 +210,7 @@ using namespace SCIRun::Core::Algorithms::Inverse;
220210

221211
// DEFINE M1 = (A * (R^T*R)^-1 * A^T MATRIX FOR FASTER COMPUTATION
222212
DenseMatrix RAtr = iRRtr * forward_transpose;
223-
M1 = *forwardMatrix_ * RAtr;
213+
M1 = forwardMatrix_ * RAtr;
224214

225215
// DEFINE M2 = (C^TC)^-1
226216
M2 = iCCtr;
@@ -232,7 +222,7 @@ using namespace SCIRun::Core::Algorithms::Inverse;
232222
M4 = DenseMatrix::Identity(M, N);
233223

234224
// DEFINE measurement vector
235-
y = *measuredData_;
225+
y = measuredData_;
236226

237227

238228

@@ -263,7 +253,7 @@ using namespace SCIRun::Core::Algorithms::Inverse;
263253
// DEFINITIONS AND PREALOCATION OF SOURCE REGULARIZATION MATRIX 'R'
264254

265255
// if R does not exist, set as identity of size equal to N (columns of fwd matrix)
266-
if (!sourceWeighting_)
256+
if (true)//(&sourceWeighting_==NULL)
267257
{
268258
RtrR = DenseMatrix::Identity(N, N);
269259
}
@@ -272,19 +262,19 @@ using namespace SCIRun::Core::Algorithms::Inverse;
272262
// if provided the non-squared version of R
273263
if( regularizationSolutionSubcase_==solution_constrained )
274264
{
275-
RtrR = sourceWeighting_->transpose() * *sourceWeighting_;
265+
RtrR = sourceWeighting_.transpose() * sourceWeighting_;
276266
}
277267
// otherwise, if the source regularization is provided as the squared version (RR^T)
278268
else if ( regularizationSolutionSubcase_==solution_constrained_squared )
279269
{
280-
RtrR = *sourceWeighting_;
270+
RtrR = sourceWeighting_;
281271
}
282272
}
283273

284274

285275
// DEFINITIONS AND PREALOCATIONS OF MEASUREMENTS COVARIANCE MATRIX 'C'
286276
// if C does not exist, set as identity of size equal to M (rows of fwd matrix)
287-
if (!sensorWeighting_)
277+
if (true)//(&sensorWeighting_==NULL)
288278
{
289279
CtrC = DenseMatrix::Identity(M, M);
290280
}
@@ -293,18 +283,18 @@ using namespace SCIRun::Core::Algorithms::Inverse;
293283
// if measurement covariance matrix provided in non-squared form
294284
if (regularizationResidualSubcase_ == residual_constrained)
295285
{
296-
CtrC = sensorWeighting_->transpose() * *sensorWeighting_;
286+
CtrC = sensorWeighting_.transpose() * sensorWeighting_;
297287
}
298288
// otherwise if the source covariance matrix is provided in squared form
299289
else if ( regularizationResidualSubcase_ == residual_constrained_squared )
300290
{
301-
CtrC = *sensorWeighting_;
291+
CtrC = sensorWeighting_;
302292
}
303293

304294
}
305295

306296
// DEFINE M1 = (A * (R*R^T)^-1 * A^T MATRIX FOR FASTER COMPUTATION
307-
DenseMatrix CtrCA = CtrC * (*forwardMatrix_);
297+
DenseMatrix CtrCA = CtrC * (forwardMatrix_);
308298
M1 = forward_transpose * CtrCA;
309299

310300
// DEFINE M2 = (CC^T)^-1
@@ -317,7 +307,7 @@ using namespace SCIRun::Core::Algorithms::Inverse;
317307
M4 = CtrCA.transpose();
318308

319309
// DEFINE measurement vector
320-
y = CtrCA.transpose() * *measuredData_;
310+
y = CtrCA.transpose() * measuredData_;
321311

322312
}
323313

‎src/Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovImpl_child.h‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -54,17 +54,19 @@ namespace SCIRun {
5454
public:
5555
SolveInverseProblemWithTikhonovImpl_child() {};
5656

57+
5758
private:
5859

5960
SCIRun::Core::Datatypes::DenseMatrix M1;
6061
SCIRun::Core::Datatypes::DenseMatrix M2;
6162
SCIRun::Core::Datatypes::DenseMatrix M3;
6263
SCIRun::Core::Datatypes::DenseMatrix M4;
63-
SCIRun::Core::Datatypes::DenseColumnMatrix y;
64+
SCIRun::Core::Datatypes::DenseMatrix y;
6465

66+
protected:
6567

66-
SCIRun::Core::Datatypes::DenseColumnMatrix computeInverseSolution( double lambda_sq, bool inverseCalculation);
67-
void preAlocateInverseMatrices(SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_,SCIRun::Core::Datatypes::DenseMatrix& measuredData_ ,SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_,SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_);
68+
virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double lambda_sq, bool inverseCalculation) const;
69+
void preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_);
6870
// bool checkInputMatrixSizes(); // DEFINED IN PARENT, MIGHT WANT TO OVERRIDE SOME OTHER TIME
6971

7072
};

‎src/Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovSVD_impl.cc‎

Lines changed: 11 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -63,15 +63,16 @@ using namespace SCIRun::Core::Algorithms;
6363
using namespace SCIRun::Core::Algorithms::Inverse;
6464

6565

66+
6667
///////////////////////////////////////////////////////////////////
6768
/////// prealocate Matrices for inverse compuation
6869
/// This function precalcualtes the SVD of the forward matrix and prepares singular vectors and values for posterior computations
6970
///////////////////////////////////////////////////////////////////
70-
void SolveInverseProblemWithTikhonovSVD_impl::preAlocateInverseMatrices( SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_,SCIRun::Core::Datatypes::DenseMatrix& measuredData_ ,SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_,SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_)
71+
void SolveInverseProblemWithTikhonovSVD_impl::preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_)
7172
{
7273

7374
// Compute the SVD of the forward matrix
74-
Eigen::JacobiSVD<SCIRun::Core::Datatypes::DenseMatrix::EigenBase> svd( *forwardMatrix_, Eigen::ComputeFullU | Eigen::ComputeFullV);
75+
Eigen::JacobiSVD<SCIRun::Core::Datatypes::DenseMatrix::EigenBase> svd( forwardMatrix_, Eigen::ComputeFullU | Eigen::ComputeFullV);
7576
SVDdecomposition = svd;
7677

7778
// determine rank
@@ -85,8 +86,7 @@ void SolveInverseProblemWithTikhonovSVD_impl::preAlocateInverseMatrices( SCIRun:
8586

8687

8788
// Compute the projection of data y on the left singular vectors
88-
SCIRun::Core::Datatypes::DenseColumnMatrix tempUy(rank);
89-
tempUy = SVDdecomposition.matrixU().transpose() * (*measuredData_);
89+
auto tempUy = SVDdecomposition.matrixU().transpose() * (measuredData_);
9090

9191
Uy = tempUy;
9292

@@ -95,22 +95,24 @@ void SolveInverseProblemWithTikhonovSVD_impl::preAlocateInverseMatrices( SCIRun:
9595
//////////////////////////////////////////////////////////////////////
9696
// THIS FUNCTION returns regularized solution by tikhonov method
9797
//////////////////////////////////////////////////////////////////////
98-
SCIRun::Core::Datatypes::DenseColumnMatrix SolveInverseProblemWithTikhonovSVD_impl::computeInverseSolution( double lambda_sq, bool inverseCalculation )
98+
SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovSVD_impl::computeInverseSolution( double lambda_sq, bool inverseCalculation ) const
9999
{
100100

101101
// prealocate matrices
102-
const int N = forwardMatrix_->ncols();
103-
const int M = forwardMatrix_->nrows();
104-
DenseColumnMatrix solution(DenseMatrix::Zero(N,1));
102+
const int N = SVDdecomposition.matrixV().cols();
103+
const int M = SVDdecomposition.matrixU().rows();
104+
const int numTimeSamples = Uy.ncols();
105+
DenseMatrix solution(DenseMatrix::Zero(N,numTimeSamples));
105106
DenseMatrix tempInverse(DenseMatrix::Zero(N,M));
106107

108+
int x = 0;
107109

108110
// Compute inverse solution
109111
for (int rr=0; rr<rank ; rr++)
110112
{
111113
// evaluate filter factor
112114
double singVal = SVDdecomposition.singularValues()[rr];
113-
double filterFactor_i = singVal / ( lambda_sq + singVal * singVal ) * Uy[rr];
115+
double filterFactor_i = singVal / ( lambda_sq + singVal * singVal ) * Uy(rr);
114116

115117
// u[date solution
116118
solution += filterFactor_i * SVDdecomposition.matrixV().col(rr);

‎src/Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovSVD_impl.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -72,11 +72,11 @@ namespace Inverse {
7272
// const SCIRun::Core::Datatypes::DenseColumnMatrix matrixS_;
7373
// const SCIRun::Core::Datatypes::DenseMatrix matrixV_;
7474
// SCIRun::Core::Datatypes::DenseMatrix y;
75-
SCIRun::Core::Datatypes::DenseColumnMatrix Uy;
75+
SCIRun::Core::Datatypes::DenseMatrix Uy;
7676

7777

78-
SCIRun::Core::Datatypes::DenseColumnMatrix computeInverseSolution( double lambda_sq, bool inverseCalculation);
79-
void preAlocateInverseMatrices(SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_,SCIRun::Core::Datatypes::DenseMatrix& measuredData_ ,SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_,SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_);
78+
virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double lambda_sq, bool inverseCalculation) const;
79+
void preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_);
8080
// bool checkInputMatrixSizes(); // DEFINED IN PARENT, MIGHT WANT TO OVERRIDE SOME OTHER TIME
8181

8282

0 commit comments

Comments
 (0)