SOFUS (SIMD Optimized Fast Ultrasound Simulator)

Author

Jens Munk Hansen

Date

2017

Table of contents

• Introduction
• Release Notes
• Build and Run
• Building Blocks
  • Mathematical primitives
  • Transducer representation
  • Frequency domain responses
• Graphical User Interface
  • Usage
• Article
• Other Stuff
  • Functional requirements
• Examples of usage
  • C++ example
  • ANSI-C example

Introduction

SOFUS is a framework and an application for simulating ultrasonic fields in frequency or time domain. The first part of this work was published in 2016 [7] . The work is made primarily is C++ and is wrapped for use in Python with the extensible compiler, simplified wrapper and interface generator (SWIG) [2] . The program is implemented to support modern single instruction, multiple data (SIMD) architectures - hence it's name SIMD Optimized Fast Ultrasound Simulator (SOFUS).

Graphical User Interface

An initial user interface has been created with Python using PyQt4. The intention is to make a native application, which is easier to distribute. For now, this is the only graphical interface to the library.

The GUI consists of a number of tabs, where each tab is used for a given purpose.

• Simulation

  This tab is used for setting simulation parameters, e.g. integration order for time or frequency domain simulations. Global parameters such as speed-of-sound is also set here.

• Transducer

  This tab is used for defining the transducer by specifying a number of parameters and using one of the supported geometries

  

• Field

  The transmitted field can be computed on any axis-parallel 2D grid. If two singleton dimensions are selected, a 1D curve will be computed.

  

Usage

Bla bla bla

Other Stuff

Functional requirements

Examples of usage

C++ example

Below, the continuous wave field is computed for a 128-element linear array.

  const float f0 = 1.0e6f;
  const float c  = 1500;
  const size_t nElements = 128;
  const float kerf   = 5.0e-4f;
  const float width  = 3e-3f;
  const float height = 50e-3f;
 
  size_t nDiv = 4;
 
  const size_t nx = 170;
  const size_t nz = 250;
 
  auto pos = sps::unique_aligned_array_create<float>(nx*nz*3);
 
  const float d = (width+kerf)*nx;
 
  const float dx = (1.5f * d) / nx;
  const float dz = (2.0f * d) / nz;
 
  float focus[3] = {0, 0, d};
 
  std::vector<float> xs(nx);
  std::vector<float> zs(nz);
 
  float wx = static_cast<float>(nx-1) / 2;
  for (size_t i = 0 ; i < nx ; i++) {
    xs[i] = (static_cast<float>(i) - wx) * dx;
  }
  for (size_t i = 0 ; i < nz ; i++) {
    zs[i] = static_cast<float>(i) * dz;
  }
 
  fnm::Aperture<float> a(nElements, width, kerf, height);
 
  a.F0Set(f0);
  a.CSet(c);
  a.NDivWSet(nDiv);
  a.NDivHSet(nDiv);
  a.NThreadsSet(4);
  a.FocusSet(focus);
  a.FocusingTypeSet(fnm::FocusingType::Rayleigh);
 
  for (size_t i = 0 ; i < nx ; i++) {
    for (size_t j = 0 ; j < nz ; j++) {
      pos[3*(i*nz+j) + 0] = xs[i];
      pos[3*(i*nz+j) + 1] = 0.0f;
      pos[3*(i*nz+j) + 2] = zs[j];
    }
  }
 
  // Output variables
  size_t nresults = 0;
  std::complex<float>* results = NULL;
 
  // Compute pressure
  int err;
  err = a.CalcCwFast(pos.get(), nx*nz, 3, &results, &nresults);
  SPS_UNREFERENCED_PARAMETER(err);
 
  // Clean-up
  free(results);

ANSI-C example

Below, the continuous wave field is computed for a 128-element linear array.

  const float f0 = 1.0e6f;
  const float c  = 1500.0f;
  const size_t nElements = 64; // 128
  const float kerf   = 5.0e-4f;
  const float width  = 3e-3f;
  const float height = 50e-3f;
 
  size_t nDiv = 4;
 
  const size_t nx = 170;
  const size_t nz = 250;
 
  float* pos = (float*) malloc(nx*nz*3*sizeof(float));
 
  const float d = (width+kerf)*nx;
 
  const float dx = (1.5f * d) / nx;
  const float dz = (2.0f * d) / nz;
 
  float focus[3] = {0,0,d};
 
  float* xs = (float*) malloc(nx*sizeof(float));
  float* zs = (float*) malloc(nz*sizeof(float));
 
  float wx = (nx-1.0f) / 2;
 
  size_t i,j;
  for (i = 0 ; i < nx ; i++) {
    xs[i] = ( (float)i - wx) * dx;
  }
  for (i = 0 ; i < nz ; i++) {
    zs[i] = (float)i * dz;
  }
 
  Aperture* a = NULL;
  int err = ApertureLinearCreate(&a, nElements, width, kerf, height);
 
#ifdef FNM_CLOSURE_FUNCTIONS
  ApertureRwFloatParamSet(a,F0,&f0,0);
  ApertureRwFloatParamSet(a,C,&c,0);
  ApertureRwFloatParamSet(a,Focus, &focus[0], 1);
#endif
  // TODO: Narrow down interface
  ApertureNDivWSet(a,nDiv);
  ApertureNDivHSet(a,nDiv);
  ApertureNThreadsSet(a,1);
 
  ApertureFocusingTypeSet(a,Rayleigh);
 
  for (i = 0 ; i < nx ; i++) {
    for (j = 0 ; j < nz ; j++) {
      pos[3*(i*nz+j) + 0] = xs[i];
      pos[3*(i*nz+j) + 1] = 0.0f;
      pos[3*(i*nz+j) + 2] = zs[j];
    }
  }
 
  size_t nresults = 0;
  float complex* results = NULL;
 
#ifndef FNM_DOUBLE_SUPPORT
  err = ApertureCalcCwFieldRef(a, pos, nx*nz, 3, &results, &nresults);
#else
  err = ApertureCalcCwFast(a, pos, nx*nz, 3, &results, &nresults);
#endif
 
  // Clean-up
  err = ApertureDestroy(a);
 
  if (results) {
    free(results);
  }
 
  free(pos);
  free(xs);
  free(zs);