utils.h

/***************************************************************************
                          mri_utils.h  -  description
                             -------------------
    begin                : Wed Feb 9 2005
    copyright            : (C) 2005 by Thies H. Jochimsen
    email                : 
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
#ifndef UTILS_H
#define UTILS_H

#include <odindata/data.h>




template <int N_rank>   
Array<float,N_rank> truncate_float(Array<float,N_rank> a){
  Array<float,N_rank> result(a.shape());
  result=where(a <= 0, ceil(a),floor(a));
  return result;
}

template <int N_rank>
void wrapPhase(Array<float,N_rank>  ph){
  ph=ph-(truncate_float(Array<float,N_rank>(ph/(2.*PII))))*2.*PII;
  ph=ph-(truncate_float(Array<float,N_rank>(ph/PII)))*2.*PII;
}

template <int N_rank>   
void unwrapPhase1d(Data<float,N_rank> ph){

    unsigned long   i,j,dim;
    float add1;
    dim=ph.extent(ph.dimensions()-1);
  Array<float,1> phu(dim);
    wrapPhase(ph);
    
  for(i=0;i<(unsigned)ph.size()/dim;i++){
    phu(0)=ph(ph.create_index(i*dim));
    add1=0;
    for(j=1;j<dim;j++){
        if ((ph(ph.create_index(i*dim+j))-ph(ph.create_index(i*dim+j-1))) > PII) add1=add1-2.*PII;
        if ((ph(ph.create_index(i*dim+j))-ph(ph.create_index(i*dim+j-1))) < -PII) add1=add1+2.*PII;
        phu(j)=ph(ph.create_index(i*dim+j))+add1;
    }
    add1=phu(dim/2);
    add1=(int)(add1/2./PII)*2.*PII+(int)((add1-(int)(add1/2./PII)*2.*PII)/PII)*2.*PII;
    for(j=0;j<dim;j++)
        ph(ph.create_index(i*dim+j))=phu(j)-add1;

        
  }
}




template<typename T>
Array<T,1> matrix_product(const Array<T,2>& matrix, const Array<T,1>& vector) {
  Log<OdinData> odinlog("","matrix_product");
  int nrows=matrix.extent()(0);
  int ncols=matrix.extent()(1);

  Array<T,1> result(nrows);
  result=0;

  int vector_extent=vector.extent(0);
  if(vector.extent(0)!=ncols) {
    ODINLOG(odinlog,errorLog) << "size mismatch (vector_extent=" << vector_extent << ") != (ncols=" << ncols << ")" << STD_endl;
    return result;
  }


  for(int icol=0; icol<ncols; icol++) {
    for(int irow=0; irow<nrows; irow++) {
      result(irow)+=matrix(irow,icol)*vector(icol);
    }
  }

  return result;
}


template<typename T>
Array<T,2> matrix_product(const Array<T,2>& matrix1, const Array<T,2>& matrix2) {
  Log<OdinData> odinlog("","matrix_product");
  int nrows=matrix1.extent()(0);
  int ncols=matrix2.extent()(1);
  ODINLOG(odinlog,normalDebug) << "nrows/ncols=" << nrows << "/" << ncols << STD_endl;

  Array<T,2> result(nrows,ncols);
  result=0;

  if(matrix1.extent(1)!=matrix2.extent(0)) {
    ODINLOG(odinlog,errorLog) << "size mismatch (matrix1=" << matrix1.shape() << ", matrix2=" << matrix2.shape() << ")" << STD_endl;
    return result;
  }

  int nprod=matrix1.extent(1);
  ODINLOG(odinlog,normalDebug) << "nprod=" << nprod << STD_endl;
  for(int irow=0; irow<nrows; irow++) {
    for(int icol=0; icol<ncols; icol++) {

      T scalprod(0);
      for(int iprod=0; iprod<nprod; iprod++) {
        scalprod+= matrix1(irow,iprod)*matrix2(iprod,icol);
      }
      result(irow,icol)=scalprod;
    }
  }

  return result;
}



template<typename T>
Array<T,1> vector_product(const Array<T,1>& u, const Array<T,1>& v) {
  Log<OdinData> odinlog("","vector_product");
  Array<T,1> result(3);
  if(u.extent(0)!=3 || v.extent(0)!=3) {
    ODINLOG(odinlog,errorLog) << "input size != 3" << STD_endl;
    return result;
  }
  result(0)=u(1)*v(2)-u(2)*v(1);
  result(1)=u(2)*v(0)-u(0)*v(2);
  result(2)=u(0)*v(1)-u(1)*v(0);
  return result;
}



template<typename T, int N_rank>
bool operator == (const TinyVector<T,N_rank>& t1, const TinyVector<T,N_rank>& t2) {
  return sum(abs(t1-t2))==0;
}


template<typename T, int N_rank>
bool operator != (const TinyVector<T,N_rank>& t1, const TinyVector<T,N_rank>& t2) {
  return !(t1==t2);
}


template<typename T, int N_rank>
bool same_shape(const Array<T,N_rank>& a1, const Array<T,N_rank>& a2) {
  return a1.shape()==a2.shape();
}




template <typename T>
bool check_range(T& val, T min, T max) {
  bool in_range=true;
  if(val<min) {val=min; in_range=false;}
  if(val>max) {val=max; in_range=false;}
  return in_range;
}




#endif

---