Adapted virtual sources to dasmtx 2d

Author: gbernardino

examples/pfield.ipynb

@@ -529,7 +529,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "gbernardino",
+   "display_name": "pymust",
    "language": "python",
    "name": "python3"
   },

examples/quickstart_demo.ipynb

@@ -1654,7 +1654,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "gbernardino",
+   "display_name": "pymust",
    "language": "python",
    "name": "python3"
   },

examples/rotatingDiskVelocitySynthetic.ipynb

@@ -7,7 +7,7 @@
    "outputs": [],
    "source": [
     "import pymust, pymust.utils\n",
-    "import numpy as np, matplotlib.pyplot as plt, scipy, scipy.io\n"
+    "import numpy as np, matplotlib.pyplot as plt, scipy, scipy.io"
    ]
   },
   {
@@ -58,13 +58,6 @@
    "execution_count": 4,
    "metadata": {},
    "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "metadata": {},
-   "outputs": [],
    "source": [
     "#% Create a 2.5-cm-by-2.5-cm image grid.\n",
     "dx = 1e-4; # grid x-step (in m)\n",
@@ -74,70 +67,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "(334, 128, 32)"
-      ]
-     },
-     "execution_count": 18,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "RF.shape"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 20,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "dtype('complex128')"
-      ]
-     },
-     "execution_count": 20,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "M.dtype"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 26,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/Users/gbernardino/pymust/src/pymust/iq2doppler.py:109: RuntimeWarning: overflow encountered in scalar multiply\n",
-      "  VN = c*PRF/4/fc/lag; #% Nyquist velocity\n",
-      "/Users/gbernardino/pymust/src/pymust/iq2doppler.py:135: RuntimeWarning: overflow encountered in scalar multiply\n",
-      "  return  param.c*PRF/4/param.fc/lag\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "682 ms ± 29 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
-    "%%timeit\n",
-    "\n",
     "# Demodulate the RF signals with RF2IQ.\n",
     "IQ =pymust.rf2iq(RF,param); \n",
     "\n",
@@ -156,42 +89,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "array([[ 1043.75942521 +393.28573606j,  -110.59674674 -794.7907838j ,\n",
-       "          155.97024474  +19.16469254j, ...,\n",
-       "          640.05718968  +79.877806j  ,    47.97989278+1650.34056168j,\n",
-       "         -786.16390502+1323.92393539j],\n",
-       "       [ 2120.85562304+1052.11123787j,   449.89643897 +473.17923728j,\n",
-       "        -2608.15759155  -45.45656041j, ...,\n",
-       "          269.26341772-1415.41186441j,  -448.39690033-2219.91092253j,\n",
-       "          268.33542192-2998.85409566j],\n",
-       "       [ -709.59873315 -475.67716058j,  -279.17991451-2085.20116294j,\n",
-       "         1562.51252753 +590.63287135j, ...,\n",
-       "        -1943.1324732  +197.70210888j,  -459.72749632 +978.46788931j,\n",
-       "         -769.83180953+1284.49579019j],\n",
-       "       ...,\n",
-       "       [ 2038.61650079-1722.36609223j,   762.66729482-1193.9384699j ,\n",
-       "        -1016.10629388 +895.42186554j, ...,\n",
-       "         2511.86589386+1001.69833543j,  1270.01473896 +886.98610086j,\n",
-       "         1889.04382837+1185.07157389j],\n",
-       "       [-3954.04195746-2266.62318821j,  -101.54895284+1101.5181209j ,\n",
-       "         1827.17382357 -666.99963509j, ...,\n",
-       "         1041.77512085-2596.96349652j,  -911.84412261-2286.95862888j,\n",
-       "        -1901.66374869-1338.63133119j],\n",
-       "       [  548.99674494+1990.93245908j,   530.02085241-1093.80990497j,\n",
-       "         -554.44475243 -914.16567814j, ...,\n",
-       "         -286.04014307-2006.75468774j,  1242.96603169 +272.22541231j,\n",
-       "         -993.40683334 +383.0611551j ]])"
-      ]
-     },
-     "execution_count": 25,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "NameError",
+     "evalue": "name 'M' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[1], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mM\u001b[49m \u001b[38;5;241m@\u001b[39m IQ\u001b[38;5;241m.\u001b[39mreshape(\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m, RF\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m2\u001b[39m])\n",
+      "\u001b[0;31mNameError\u001b[0m: name 'M' is not defined"
+     ]
     }
    ],
    "source": [

examples/rotatingDiskVelocitySyntheticML.ipynb

@@ -1,4 +1,6 @@
 import numpy as np
+import scipy.optimize, itertools 
+
 import scipy, scipy.interpolate
 from . import  utils
 
@@ -439,10 +441,56 @@ def dasmtx(SIG : np.ndarray, x: np.ndarray, z: np.ndarray, *varargin):
     #% For a given location [x(k),z(k)], one has:
     #% dTX(k) = min(delaysTXi*c + sqrt((xTi-x(k)).^2 + (zTi-z(k)).^2));
     #% In a compact matrix form, this yields:
-    if not param.passive:
-        dTX = np.min(delaysTXi*c + np.sqrt((xTi-x)**2 + (zTi-z)**2), 1).reshape((-1,1))
-    else:
+    if param.passive:
         dTX = np.zeros_like(x)
+    else:
+        # OLD version of computing the transmission delays
+        #% For a given location [x(k),z(k)], one has:
+        #% dTX(k) = min(delaysTXi*c + sqrt((xTi-x(k)).^2 + (zTi-z(k)).^2));
+        #% In a compact matrix form, this yields:
+        #dTX = np.min(delaysTXi*c + np.sqrt((xTi-x)**2 + (zTi-z)**2), 1).reshape((-1,1))
+
+        WasTransmitting = np.logical_not(np.isnan(delaysTX))
+        assert np.sum(np.abs(np.diff(WasTransmitting)))<3, 'Multiple transmitting sub-apertures are not allowed during beamforming.'
+        nTX = np.count_nonzero(WasTransmitting); # number of transmitting elements
+
+
+        #-- TX distances
+        if nTX==1:
+            # Only one element was active
+            dTX = np.hypot(xe[WasTransmitting][0]-x,ze[WasTransmitting][0]-z) + delaysTX[WasTransmitting][0]*c
+
+        elif nTX<3:
+            raise ValueError('If not 1, the number of neighboring transmitting elements must be at least 3.')
+        
+        else:
+            #-- We create a virtual transducer to calculate the TX distances.
+            xv,zv = vxdcr(xe[WasTransmitting],ze[WasTransmitting], delaysTX[WasTransmitting],c);
+            dzv = diff2(xv,zv)
+
+            # Distances between the (x,z) points and the lines normal to the
+            # virtual transducer at (xv,zv)
+            Dn = np.abs(x-xv + dzv*(z-zv))/np.hypot(1,dzv);
+            # idx contains the element numbers that give the smallest Dn_s
+            idx = np.argmin(Dn,1)
+            
+            InterpMethod = 'nearest';
+
+            if InterpMethod == 'nearest':
+                dTX = np.hypot(xv[idx].reshape((-1,1)) - x, zv[idx].reshape((-1,1))  -z)
+            elif InterpMethod == 'parabolic':
+                N = x.size
+                dTX = np.hypot(xv.reshape((1,-1))-x,zv.reshape((1,-1))-z);
+                dTX[:,1:nc] = scipy.signal.convolve2d(dTX, np.array([1, 0, 1]).reshape((1,-1)),'valid')/2 - \
+                    scipy.signal.convolve2d(dTX,np.array([1, 0, -1]).reshape((1,-1)),'valid')* \
+                    scipy.signal.convolve2d(Dn,np.array([1, 0, -1]).reshape((1,-1)),'valid')/ \
+                scipy.signal.convolve2d(Dn,np.array([2, -4, 2]).reshape((1,-1)),'valid')
+
+                dTX = dTX.flatten()
+                dTX = dTX[np.arange(N) + idx *N ]
+            else:
+                raise ValueError("Incorrect interpolation method on dTX. Valid options: [nearest | parabolic]")
+
 
 
     #%-- RX distances
@@ -588,4 +636,65 @@ def __apply__(self, SIG, axis ):
         return (self.M @ SIG.flatten(order = 'F')).reshape(self.shape, order = 'F')
 
 
-    
+# Function to calculate the first derivative of y with respect to x using a second-order difference method
+def diff2(x, y):
+    m = len(y)
+    dx = np.diff(x)
+    dy = np.zeros_like(y)
+
+    # y'(x_1)
+    dy[0] = (1 / dx[0] + 1 / dx[1]) * (y[1] - y[0]) + dx[0] / (dx[0] * dx[1] + dx[1]**2) * (y[0] - y[2])
+
+    # y'(x_m)
+    dy[-1] = (1 / dx[-2] + 1 / dx[-1]) * (y[-1] - y[-2]) + dx[-1] / (dx[-1] * dx[-2] + dx[-2]**2) * (y[-3] - y[-1])
+
+    # y'(x_i) (i > 1 & i < m)
+    dx1 = dx[:-1]
+    dx2 = dx[1:]
+    y1 = y[:-2]
+    y2 = y[1:-1]
+    y3 = y[2:]
+    dy[1:-1] = 1 / (dx1 * dx2 * (dx1 + dx2)) * (-dx2**2 * y1 + (dx2**2 - dx1**2) * y2 + dx1**2 * y3)
+
+    return dy
+
+def vxdcr(xe, ze, delaysTX, c):
+    # Virtual transducer (XDCR)
+    T2 = delaysTX**2
+    dT2dxe = diff2(xe, T2)
+    err = np.min(np.diff(xe)) * 1e-3
+    dzedxe = diff2(xe, ze)
+
+    # Virtual transducer positions
+    if np.all(ze == 0):
+        # For a linear transducer
+        xV = xe - 0.5 * c**2 * dT2dxe
+        zV = -np.sqrt(np.abs(c**2 * T2 - (xV - xe)**2))
+    else:
+        # For a convex array
+        xV = xe.copy()
+        tmp = np.sqrt(np.abs(c**2 * T2 - (xV - xe)**2))
+
+        def myfun(xv, k):
+            return (xe[k] - xv) + tmp[k] * dzedxe[k] - 0.5 * c**2 * dT2dxe[k]
+
+        for k in range(len(xV)):
+            xV[k] = scipy.optimize.fsolve(myfun, xV[k], args=(k,), xtol=err)[0]
+        zV = ze - tmp
+
+    # Validate the virtual transducer
+    test = np.all(np.diff(xV) > 0) and np.all((c**2 * T2 - (xV - xe)**2) > -err)
+    if not test:
+        raise ValueError("The delays do not allow a virtual transducer to be calculated.")
+
+    # Check if the virtual transducer is convex or concave
+    idx = np.array([list(i) for i in itertools.combinations(range(len(xV)), 3)])
+    lambda_vals = (xV[idx[:, 2]] - xV[idx[:, 1]]) / (xV[idx[:, 0]] - xV[idx[:, 1]])
+    alst1 = np.abs(lambda_vals) < 1
+    tmp = zV[idx[alst1, 2]] - (lambda_vals[alst1] * zV[idx[alst1, 0]] + (1 - lambda_vals[alst1]) * zV[idx[alst1, 1]])
+    test = max(np.count_nonzero(tmp > -np.finfo(float).eps),
+               np.count_nonzero(tmp < np.finfo(float).eps)) / len(tmp) > 0.999
+    if not test:
+        print("Warning: The virtual transducer is not convex or concave.")
+
+    return xV, zV

examples/rotatingDisk_real.ipynb

@@ -863,7 +863,7 @@ def pfield(x : np.ndarray,y : np.ndarray, z: np.ndarray, delaysTX : np.ndarray,
 
     #% RMS acoustic pressure (if we are in PFIELD only)
     if not (isSIMUS or isMKMOVIE):
-        RP =np.sqrt(RP).reshape(siz0)
+        RP =np.sqrt(RP).reshape(siz0, order = 'F')
         SPECT = np.swapaxes(SPECT, 0, 1)
         SPECT = SPECT.reshape([siz0[0], siz0[1], nSampling], order = 'F')
     return RP, SPECT, IDX

src/pymust/dasmtx.py

@@ -72,6 +72,15 @@
     "If this does not work, you can copy the source code from the github repository, and install it locally."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "!pip install pymust"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -652,11 +661,9 @@
     "xScatterers = np.array([0, 0,np.sqrt(2)*1e-2, - np.sqrt(2)*1e-2,  0])\n",
     "zScatterers = np.array([1e-2, 2e-2, np.sqrt(2) * 1e-2, np.sqrt(2) * 1e-2, 2.2e-2])\n",
     "RCScatterers = np.random.rand(5) # There are 5 scatterers\n",
-    "#RCScatterers[1:] = 0\n",
     "activationDelaysTX = np.array([0]).reshape(1,1)\n",
     "\n",
     "RF, _ = pymust.simus(xScatterers, zScatterers, RCScatterers, activationDelaysTX, param_single_element)\n",
-    "RF, c = pymust.tgc(RF)\n",
     "t = np.arange(RF.shape[0]) / param_single_element.fs # param_single_element.fs is the sampling frequency, needed to go from sample number to time\n"
    ]
   },