Biotechnology Industry

Molecular Imaging

ART (Advanced Research Technologies) Inc., develops and commercializes molecular imaging products for the healthcare and pharmaceutical industries.

Project Type: Mathematical programming (physics modeling)

Refactored a package with a modular approach (fit and optimization problem), Optimization (the algorithm was more than twice as fast as the original).

Client: ART Inc.

Technologies: C++, OO, Matlab, gnuplot, VTK visualizing tool, numerical programming, non-linear fit and optimization

Role: Scientific Programmer

I was hired as a physics programmer for programming physics algorithm. I have re-designed the curve fit library (optimization problem). Particularly, I have implemented an Object-Oriented architecture, modular approach by developing an easy-to-use programming API. This new version of the library featured: efficient and robust optimization for fast floating point operation, extensible class hierarchies that let user readily implement new algorithms.

Molecular Imaging Technology

Molecular imaging is a technology based on propagation of light (propagating in a medium, human body), by analyzing the optical property of light (near infrared), to diagnostic cancer.

Molecular imaging Technology is based on the diffusion of light in human body (Light propagate through human body).

This enables physicians to peer into the living body in order to identify diseases, monitor their progression or treat medical conditions at a molecular level. 

Molecular imaging involves illuminating the part of the body of interest with, non-ionizing radiation and analyzing the emergent pattern for signs of pathology. Task consist of investigating both the scattering and absorption properties of tissue, the manner in which light propagates through tissue (normal tissue differs from abnormal in its absorption or scattering), it then becomes possible to differentiate optically between normal and abnormal conditions. Molecular imaging techniques are well suited for the evaluation of the response to therapy in cancer.

Need to solve the diffusion equation by some approximating methods, called numerical method (exact solution is not known we approximate by a solution that is close the real one according to some assumptions). In this case we use an optimization method which consist (numerical approx. of the equation (non-linear fit data to find parameters))

algorithm development and implementation and programming support of scientific research.

Architecture API

Present some UML (Unified Modeling Language) of the API of the numerical library. UML diagrams are blueprint of the architecture, divide in mainly in 3 types: Use case, static class and sequence. The former … at the beginning of the project when defining requirements (grasp how its going to be use, major impact on the class design). Once this is done, ready layout classes and dependency. Finally sequence diagram show the trace … to be completed interaction between components.

Code Snippets

Top level class (abstract) which means we cannot instantiate. Provides basic services to solve an optimization problem (by some optimization method which is an attribute). The optimization method is pass as an argument to constructor (user have to create a method).

//Abstract Base Class for Curves Fitting Solver 
class OptimizationSolver {

  public:
      // Default ctor
      OptimizationSolver( std::vector<float> aVec, uint32 aMaxIter);
      OptimizationSolver( std::vector<float> aAccuracy, /*eps*/ uint32 aMaxIter, 
                          boost::shared_ptr&ltOptimizationMethod&gt aOpMeth);

      // Force destructor to be virtual
      virtual ~OptimizationSolver() {/*nothing to do*/}

      // Solve methods
      void solve( ObjectiveFunction & aFunc, std::vector<float> aGuess, std::vector<float> aStep) const;
      void solve( ObjectiveFunction & aFunc, std::vector<float> aGuess, std::vector<float> aMin, 
                  std::vector<float> aMax) const;

      // Boundary of the parameters
      void setLowBound( std::vector<float> aLowBound) { mLowBound = aLowBound;}
      void setHighBound( std::vector<float> aHighBound) { mHiBound = aHighBound;}
      std::vector<float> getLowBound() const { return mLowBound;}
      std::vector<float> getHighBound() const { return mHiBound;}

      // Initialize parameters .... and overwrited in the subclass
      virtual bool initialize( std::vector<float> aGuess, std::vector<float> aStep);

  protected:
      // Initial guess
      mutable std::vector<float> mRoot;
      mutable std::vector<float> mStep;

      // Required accuracy of the solver
      std::vector<float> mAccuracy;  // bestAccuracy_;
  
      // Maximum and real number of iterations
      unsigned int mMaxIterations;
      unsigned int mNbIterations;

      // Maximum and evaluation number of ....
      unsigned int mMaxEvaluations;
      mutable unsigned int mEvaluationNumber;

      // This method must be implemented in derived classes and contains the actual code which 
      // searches for the zeroes of the Objective Function. It assumes that:
      //  - xMin and xMax form a valid bracket
      //  - ... to be completed
      //  - mRoot was initialized to a valid initial guess
      virtual void solve( ObjectiveFunction & aFunc, std::vector<float> aAccuracy) const = 0;

      // Optimization method
      boost::shared_ptr&ltOptimizationMethod&gt mOptMeth;

  private:
      // Boundary value  
      std::vector<float> mLowBound;
      std::vector<float> mHiBound;

      // Boundary parameters almost fixed
      bool mLowBoundEnforced;
      bool mHiBoundEnforced;

      // Fix the lower and upper bounds of ...
      void enforcebounds( std::vector<float> aBounds) const;
};

Below inheriting class from the base class (depends on the fit algorithm choosen)

// "Model of" of the two-steps fitting algorithm. 
class TwoStepsFitSolver : public OptimizationSolver 
{
   public:
       // Set Direct Search algorithm as default Optimization Method
       TwoStepsFitSolver( std::vector<float> aAccuracy /*eps*/, uint32 aMaxIter, uint32 aNbParams);

       // Default ctor
       TwoStepsFitSolver( std::vector<float> aAccuracy /*eps*/, uint32 aMaxIter, uint32 aNbParams,
                          boost::shared_ptr&ltOptimizationMethod&gt aOptMeth);

       // Destructor
       ~TwoStepsFitSolver() {/*nothing to do*/}

       // Parameters of the fit
       inline ART::Data::BackGParams2Fit getParams2Fit() 
       {
           float* wResults = mOptMeth->GetSolution();
           mParams2Fit.mMua   = wResults[0];
           mParams2Fit.mMusp  = wResults[1];
           mParams2Fit.mAmpl  = wResults[2];
           mParams2Fit.mTzero = wResults[3];

           return mParams2Fit;
       }

       // Interface
       // Override solve method from base class
       void solve( ObjectiveFunction & aFunc, std::vector<float> aAccuracy) const;

       // Derived from base class
       bool initialize( std::vector<float> aGuess, std::vector<float> aStep);

   private:
       // Number of parameters for the fit
       unsigned int mNbParams2Fit;
       ART::Data::BackGParams2Fit mParams2Fit;
       // Store the result
       float* mResults;
};

Below model of a function using “called functor” in C++. More description coming

// Contini term structure function is the following  
//     T(rho, t) = exp()/2.(4piDv) Sum(......)
//
//     We evaluate the theoretical model, in this case we use the Contini Model.
//     Reference can be found in the following paper: 
//     "Photon Migration through a turbid slab described by a Model based
//      on Diffusion Approximation: I Theory", D. Contini, F. Martelli and G. Zaccanti
//     Applied Optics, Vol. 36, No.19 (Jul. 1997)
//
//  "Model of" functor
class ContiniFunction 
{
  public:
      // Default constructor to set coeffs
      ContiniFunction( float aCurrentMua, float aCurrentMusp, float aSpeedOfLight, float aThickness,
                       const ART::Data::Timing aTpsfTimeParams);
      // Operator to evaluate the function at given time
      ART::Data::FloatVec operator() ( float aRho, const float aTime);
      // Scale for ...
      void setScaleZb( float aScaleZb) { mScaledZB = aScaleZb;} 
   private:  
      // Coefficients of the function
      float mCurrentMua;
      float mCurrentMusp;
      float mSpeedOfLight;
      float mThickness;
      float mScaledZB;

      // Time parameters for the TPSF (Temporal ... Spread Function)
      ART::Data::Timing mTPSFTime;
};

Technical report “Background Fit Algorithm”, – J. Belanger, Technical Note ART Inc. (June 2005)