Skip to content

Commit

Permalink
Maint: restructure/re-organize HpiFit.fit()
Browse files Browse the repository at this point in the history
  • Loading branch information
@RDoerfel @juangpc
RDoerfel authored and juangpc committed Feb 23, 2022
1 parent dd68033 commit cf25a47
Show file tree
Hide file tree
Showing 2 changed files with 172 additions and 182 deletions.
291 changes: 143 additions & 148 deletions libraries/inverse/hpiFit/hpifit.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,11 +143,19 @@ void HPIFit::fit(const MatrixXd& matProjectedData,
const MatrixXd matAmplitudes = computeAmplitudes(matProjectedData,
hpiModelParameters);

CoilParam fittedCoilParams = computeCoilLocation(matAmplitudes,
matProjectors,
hpiFitResult.devHeadTrans,
hpiFitResult.errorDistances,
matCoilsHead);
const MatrixXd matCoilsSeed = computeSeedPoints(matAmplitudes,
hpiFitResult.devHeadTrans,
hpiFitResult.errorDistances,
matCoilsHead);

CoilParam fittedCoilParams = dipfit(matCoilsSeed,
m_sensors,
matAmplitudes,
hpiModelParameters.iNHpiCoils(),
matProjectors,
500,
1e-9f);

if(bOrderFrequencies) {
const std::vector<int> vecOrder = findCoilOrder(fittedCoilParams.pos,
matCoilsHead);
Expand Down Expand Up @@ -202,15 +210,12 @@ Eigen::MatrixXd HPIFit::computeAmplitudes(const Eigen::MatrixXd& matProjectedDat

//=============================================================================================================

CoilParam HPIFit::computeCoilLocation(const Eigen::MatrixXd& matAmplitudes,
const MatrixXd& matProjectors,
const FIFFLIB::FiffCoordTrans& transDevHead,
const QVector<double>& vecError,
const MatrixXd& matCoilsHead,
const int iMaxIterations,
const float fAbortError)
Eigen::MatrixXd HPIFit::computeSeedPoints(const Eigen::MatrixXd& matAmplitudes,
const FIFFLIB::FiffCoordTrans& transDevHead,
const QVector<double>& vecError,
const Eigen::MatrixXd& matCoilsHead)
{
const int iNumCoils = matAmplitudes.cols();
const int iNumCoils = matCoilsHead.rows();
MatrixXd matCoilsSeed = MatrixXd::Zero(iNumCoils,3);

const double dError = std::accumulate(vecError.begin(), vecError.end(), .0) / vecError.size();
Expand Down Expand Up @@ -240,94 +245,76 @@ CoilParam HPIFit::computeCoilLocation(const Eigen::MatrixXd& matAmplitudes,
}
}
}

// dipole fit
return dipfit(matCoilsSeed,
m_sensors,
matAmplitudes,
iNumCoils,
matProjectors,
iMaxIterations,
fAbortError);

return matCoilsSeed;
}

//=============================================================================================================

Eigen::VectorXd HPIFit::computeGoF(const Eigen::VectorXd& vecDipFitError)
{
VectorXd vecGoF(vecDipFitError.size());
for(int i = 0; i < vecDipFitError.size(); ++i) {
vecGoF(i) = 1 - vecDipFitError(i);
}
return vecGoF;
}

//=============================================================================================================

FIFFLIB::FiffCoordTrans HPIFit::computeDeviceHeadTransformation(const Eigen::MatrixXd& matCoilsDev,
const Eigen::MatrixXd& matCoilsHead)
CoilParam HPIFit::dipfit(const MatrixXd matCoilsSeed,
const SensorSet& sensors,
const MatrixXd& matData,
const int iNumCoils,
const MatrixXd& matProjectors,
const int iMaxIterations,
const float fAbortError)
{
MatrixXd matTrans = computeTransformation(matCoilsHead,matCoilsDev);
return FiffCoordTrans::make(1,4,matTrans.cast<float>(),true);
}

//=============================================================================================================
//Do this in conncurrent mode
//Generate QList structure which can be handled by the QConcurrent framework
QList<HPIFitData> lCoilData;

QVector<double> HPIFit::computeEstimationError(const Eigen::MatrixXd& matCoilsDev,
const Eigen::MatrixXd& matCoilsHead,
const FIFFLIB::FiffCoordTrans& transDevHead)
{
//Calculate Error
MatrixXd matTemp = matCoilsDev;
MatrixXd matTestPos = transDevHead.apply_trans(matTemp.cast<float>()).cast<double>();
MatrixXd matDiffPos = matTestPos - matCoilsHead;
for(qint32 i = 0; i < iNumCoils; ++i) {
HPIFitData coilData;
coilData.m_coilPos = matCoilsSeed.row(i);
coilData.m_sensorData = matData.col(i);
coilData.m_sensors = sensors;
coilData.m_matProjector = matProjectors;
coilData.m_iMaxIterations = iMaxIterations;
coilData.m_fAbortError = fAbortError;

// compute error
int iNumCoils = matCoilsDev.rows();
QVector<double> vecError(iNumCoils);
for(int i = 0; i < matDiffPos.rows(); ++i) {
vecError[i] = matDiffPos.row(i).norm();
lCoilData.append(coilData);
}
return vecError;
}
//Do the concurrent filtering
CoilParam coil(iNumCoils);

//=============================================================================================================
if(!lCoilData.isEmpty()) {
// //Do sequential
// for(int l = 0; l < lCoilData.size(); ++l) {
// doDipfitConcurrent(lCoilData[l]);
// }

FIFFLIB::FiffDigPointSet HPIFit::getFittedPointSet(const Eigen::MatrixXd& matCoilsDev)
{
FiffDigPointSet fittedPointSet;
int iNumCoils = matCoilsDev.rows();
//Do concurrent
QFuture<void> future = QtConcurrent::map(lCoilData,
&HPIFitData::doDipfitConcurrent);
future.waitForFinished();

for(int i = 0; i < iNumCoils; ++i) {
FiffDigPoint digPoint;
digPoint.kind = FIFFV_POINT_EEG; //Store as EEG so they have a different color
digPoint.ident = i;
digPoint.r[0] = matCoilsDev(i,0);
digPoint.r[1] = matCoilsDev(i,1);
digPoint.r[2] = matCoilsDev(i,2);

fittedPointSet << digPoint;
//Transform results to final coil information
for(qint32 i = 0; i < lCoilData.size(); ++i) {
coil.pos.row(i) = lCoilData.at(i).m_coilPos;
coil.mom = lCoilData.at(i).m_errorInfo.moment.transpose();
coil.dpfiterror(i) = lCoilData.at(i).m_errorInfo.error;
coil.dpfitnumitr(i) = lCoilData.at(i).m_errorInfo.numIterations;

//std::cout<<std::endl<< "HPIFit::dipfit - Itr steps for coil " << i << " =" <<coil.dpfitnumitr(i);
}
}
return fittedPointSet;
return coil;
}


//=============================================================================================================

std::vector<int> HPIFit::findCoilOrder(const MatrixXd& matCoilsDev,
const MatrixXd& matCoilsHead)
{
// extract digitized and fitted coils
MatrixXd matDigTemp = matCoilsHead;
MatrixXd matCoilTemp = matCoilsDev;
int iNumCoils = matCoilsDev.rows();
const int iNumCoils = matCoilsDev.rows();

std::vector<int> vecOrder(iNumCoils);
std::iota(vecOrder.begin(), vecOrder.end(), 0);;

// maximum 10 mm mean error
double dErrorMin = 0.010;
const double dErrorMin = 0.010;
double dErrorActual = 0.0;
double dErrorBest = dErrorMin;

Expand All @@ -354,10 +341,38 @@ std::vector<int> HPIFit::findCoilOrder(const MatrixXd& matCoilsDev,

//=============================================================================================================

double HPIFit::objectTrans(const MatrixXd& matHeadCoil,
const MatrixXd& matCoil,
const MatrixXd& matTrans)
{
// Compute the fiducial registration error - the lower, the better.
const int iNumCoils = matHeadCoil.rows();
MatrixXd matTemp = matCoil;

// homogeneous coordinates
matTemp.conservativeResize(matCoil.rows(),matCoil.cols()+1);
matTemp.block(0,3,iNumCoils,1).setOnes();
matTemp.transposeInPlace();

// apply transformation
MatrixXd matTestPos = matTrans * matTemp;

// remove
MatrixXd matDiff = matTestPos.block(0,0,3,iNumCoils) - matHeadCoil.transpose();
VectorXd vecError = matDiff.colwise().norm();

// compute error
double dError = matDiff.colwise().norm().mean();;

return dError;
}

//=============================================================================================================

Eigen::MatrixXd HPIFit::order(const std::vector<int>& vecOrder,
const Eigen::MatrixXd& matToOrder)
{
int iNumCoils = vecOrder.size();
const int iNumCoils = vecOrder.size();
MatrixXd matToOrderTemp = matToOrder;

for(int i = 0; i < iNumCoils; i++) {
Expand All @@ -371,7 +386,7 @@ Eigen::MatrixXd HPIFit::order(const std::vector<int>& vecOrder,
QVector<int> HPIFit::order(const std::vector<int>& vecOrder,
const QVector<int>& vecToOrder)
{
int iNumCoils = vecOrder.size();
const int iNumCoils = vecOrder.size();
QVector<int> vecToOrderTemp = vecToOrder;

for(int i = 0; i < iNumCoils; i++) {
Expand All @@ -382,55 +397,22 @@ QVector<int> HPIFit::order(const std::vector<int>& vecOrder,

//=============================================================================================================

CoilParam HPIFit::dipfit(const MatrixXd matCoilsSeed,
const SensorSet& sensors,
const MatrixXd& matData,
const int iNumCoils,
const MatrixXd& matProjectors,
const int iMaxIterations,
const float fAbortError)
Eigen::VectorXd HPIFit::computeGoF(const Eigen::VectorXd& vecDipFitError)
{
//Do this in conncurrent mode
//Generate QList structure which can be handled by the QConcurrent framework
QList<HPIFitData> lCoilData;

for(qint32 i = 0; i < iNumCoils; ++i) {
HPIFitData coilData;
coilData.m_coilPos = matCoilsSeed.row(i);
coilData.m_sensorData = matData.col(i);
coilData.m_sensors = sensors;
coilData.m_matProjector = matProjectors;
coilData.m_iMaxIterations = iMaxIterations;
coilData.m_fAbortError = fAbortError;

lCoilData.append(coilData);
VectorXd vecGoF(vecDipFitError.size());
for(int i = 0; i < vecDipFitError.size(); ++i) {
vecGoF(i) = 1 - vecDipFitError(i);
}
//Do the concurrent filtering
CoilParam coil(iNumCoils);

if(!lCoilData.isEmpty()) {
// //Do sequential
// for(int l = 0; l < lCoilData.size(); ++l) {
// doDipfitConcurrent(lCoilData[l]);
// }

//Do concurrent
QFuture<void> future = QtConcurrent::map(lCoilData,
&HPIFitData::doDipfitConcurrent);
future.waitForFinished();

return vecGoF;
}

//Transform results to final coil information
for(qint32 i = 0; i < lCoilData.size(); ++i) {
coil.pos.row(i) = lCoilData.at(i).m_coilPos;
coil.mom = lCoilData.at(i).m_errorInfo.moment.transpose();
coil.dpfiterror(i) = lCoilData.at(i).m_errorInfo.error;
coil.dpfitnumitr(i) = lCoilData.at(i).m_errorInfo.numIterations;
//=============================================================================================================

//std::cout<<std::endl<< "HPIFit::dipfit - Itr steps for coil " << i << " =" <<coil.dpfitnumitr(i);
}
}
return coil;
FIFFLIB::FiffCoordTrans HPIFit::computeDeviceHeadTransformation(const Eigen::MatrixXd& matCoilsDev,
const Eigen::MatrixXd& matCoilsHead)
{
const MatrixXd matTrans = computeTransformation(matCoilsHead,matCoilsDev);
return FiffCoordTrans::make(1,4,matTrans.cast<float>(),true);
}

//=============================================================================================================
Expand Down Expand Up @@ -503,6 +485,46 @@ Eigen::Matrix4d HPIFit::computeTransformation(Eigen::MatrixXd matNH, MatrixXd ma

//=============================================================================================================

QVector<double> HPIFit::computeEstimationError(const Eigen::MatrixXd& matCoilsDev,
const Eigen::MatrixXd& matCoilsHead,
const FIFFLIB::FiffCoordTrans& transDevHead)
{
//Calculate Error
MatrixXd matTemp = matCoilsDev;
MatrixXd matTestPos = transDevHead.apply_trans(matTemp.cast<float>()).cast<double>();
MatrixXd matDiffPos = matTestPos - matCoilsHead;

// compute error
int iNumCoils = matCoilsDev.rows();
QVector<double> vecError(iNumCoils);
for(int i = 0; i < matDiffPos.rows(); ++i) {
vecError[i] = matDiffPos.row(i).norm();
}
return vecError;
}

//=============================================================================================================

FIFFLIB::FiffDigPointSet HPIFit::getFittedPointSet(const Eigen::MatrixXd& matCoilsDev)
{
FiffDigPointSet fittedPointSet;
const int iNumCoils = matCoilsDev.rows();

for(int i = 0; i < iNumCoils; ++i) {
FiffDigPoint digPoint;
digPoint.kind = FIFFV_POINT_EEG; //Store as EEG so they have a different color
digPoint.ident = i;
digPoint.r[0] = matCoilsDev(i,0);
digPoint.r[1] = matCoilsDev(i,1);
digPoint.r[2] = matCoilsDev(i,2);

fittedPointSet << digPoint;
}
return fittedPointSet;
}

//=============================================================================================================

void HPIFit::storeHeadPosition(float fTime,
const Eigen::MatrixXf& transDevHead,
Eigen::MatrixXd& matPosition,
Expand All @@ -528,30 +550,3 @@ void HPIFit::storeHeadPosition(float fTime,
matPosition(matPosition.rows()-1,9) = 0;
}

//=============================================================================================================

double HPIFit::objectTrans(const MatrixXd& matHeadCoil,
const MatrixXd& matCoil,
const MatrixXd& matTrans)
{
// Compute the fiducial registration error - the lower, the better.
int iNumCoils = matHeadCoil.rows();
MatrixXd matTemp = matCoil;

// homogeneous coordinates
matTemp.conservativeResize(matCoil.rows(),matCoil.cols()+1);
matTemp.block(0,3,iNumCoils,1).setOnes();
matTemp.transposeInPlace();

// apply transformation
MatrixXd matTestPos = matTrans * matTemp;

// remove
MatrixXd matDiff = matTestPos.block(0,0,3,iNumCoils) - matHeadCoil.transpose();
VectorXd vecError = matDiff.colwise().norm();

// compute error
double dError = matDiff.colwise().norm().mean();;

return dError;
}
Loading

0 comments on commit cf25a47

Please sign in to comment.