utils.h

/***************************************************************************
                          mri_utils.h  -  description
                             -------------------
    begin                : Wed Feb 9 2005
    copyright            : (C) 2000-2021 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>
 
Data<float,1> unwrap_phase(const Data<float,1>& phase, int startindex=0);
 
 
 
 
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_rows, int N_columns>
TinyVector<T,N_rows> operator * (const TinyMatrix<T, N_rows, N_columns>& matrix, const TinyVector<T, N_columns>& vector) {
  TinyVector<T,N_rows> result;
  result=0;
 
  for(int icol=0; icol<N_columns; icol++) {
    for(int irow=0; irow<N_rows; irow++) {
      result(irow)+=matrix(irow,icol)*vector(icol);
    }
  }
 
  return result;
}
 
 
template<typename T, int N_rank>
bool same_shape(const Array<T,N_rank>& a1, const Array<T,N_rank>& a2, const TinyVector<int,N_rank>& dimmask=1) {
  for(int i=0; i<N_rank; i++) {
    if(dimmask(i) && (a1.extent(i)!=a2.extent(i)) ) return false;
  }
  return true;
}
 
 
template<int N_rank>
bool on_grid(const TinyVector<int,3>& shape, const TinyVector<int,3>& index) {
  for(int i=0; i<N_rank; i++) {
    if(index(i)<0 || index(i)>=shape(i)) {
      return false;
    }
  }
  return true;
}
 
 
 
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;
}
 
 
 
template<typename T, int N_rank>
void clip_max(Data<T,N_rank>& data, T val) {
  TinyVector<int,N_rank> indexvec;
  for(unsigned int i=0; i<data.numElements(); i++) {
    indexvec=data.create_index(i);
    if(data(indexvec)>val) data(indexvec)=val;
  }
}
 
template<typename T, int N_rank>
void clip_min(Data<T,N_rank>& data, T val) {
  TinyVector<int,N_rank> indexvec;
  for(unsigned int i=0; i<data.numElements(); i++) {
    indexvec=data.create_index(i);
    if(data(indexvec)<val) data(indexvec)=val;
  }
}
 
 
 
#endif