jacobian.h

// -*-C++-*-
// =========================================================================
// jacobian.h (libstoast)
// Interface for Jacobian calculation
//
// Jacobian of forward operator: derivative of measurables with respect to
// parameters (frequency domain case)
// The Jacobian is a real matrix is composed of 4 blocks:
//
//     |  dlnamp  |  dlnamp  |
//     |  ------  |  ------  |
//     |   dmua   |  dkappa  |
//     |----------+----------|
//     |  dphase  |  dphase  |
//     |  ------  |  ------  |
//     |   dmua   |  dkappa  |
//
// Each row corresponds to a measurement (source-detector pair), each column
// corresponds to a voxel in solution space.
// The returned Jacobian does not contain any data or parameter scalings.
// =========================================================================

#ifndef __JACOBIAN_H
#define __JACOBIAN_H

#include "mathlib.h"
#include "felib.h"
#include "toasttype.h"

// Flags for computing Jacobian on nodal/element basis rather than
// basis provided by a raster object
#define INT_NDBASIS 0
#define INT_ELBASIS -1
#define RASTER_NDBASIS ((Raster*)INT_NDBASIS)
#define RASTER_ELBASIS ((Raster*)INT_ELBASIS)


STOASTLIB void GenerateJacobian (const Raster *raster, const QMMesh *mesh,
    const CCompRowMatrix &mvec, const CVector *dphi, const CVector *aphi,
    DataScale dscale, RDenseMatrix &J);

STOASTLIB void GenerateJacobian (const Raster *raster, const QMMesh *mesh,
    const CVector *dphi, const CVector *aphi, const CVector *proj,
    DataScale dscale, RDenseMatrix &J);

STOASTLIB void GenerateJacobian_cw (const Raster *raster, const QMMesh *mesh,
    const RVector *dphi, const RVector *aphi,
    RDenseMatrix *Jmua = NULL, RDenseMatrix *Jkap = NULL);

STOASTLIB void GenerateJacobian_cw (const Raster *raster,
    const QMMesh *mesh, const RCompRowMatrix &mvec,
    const RVector *dphi, const RVector *aphi,
    DataScale dscale, RDenseMatrix *Jmua = NULL, RDenseMatrix *Jkap = NULL);

STOASTLIB void GenerateJacobian_cw (const Raster *raster,
    const QMMesh *mesh, const RVector *dphi, const RVector *aphi,
    const RVector *proj, DataScale dscale,
    RDenseMatrix *Jmua = NULL, RDenseMatrix *Jkap = NULL);


template<class T>
STOASTLIB void ImageGradient (const IVector &dim, const RVector &size,
    const TVector<T> &im, TVector<T> *grad, const int *mask = 0);

#endif // !__JACOBIAN_H