diff --git a/code/part4_spectral_ct.ipynb b/code/part4_spectral_ct.ipynb index 0c67d5c..bb3dbf4 100644 --- a/code/part4_spectral_ct.ipynb +++ b/code/part4_spectral_ct.ipynb @@ -80,30 +80,9 @@ }, { "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEKCAYAAADTgGjXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8nGW5+P/PlclMZpImadKmpfsCpbSF0kIocNhFSgGBoigU9RTxyEFw/+Gx6hEFRQH5uZ1yUOQg4kEKVMQeFRFZVBDpAqVAF7rShpa2Wdo0yUwymbm+fzzPTCfpJJlJMkuS6/16zSszzzb3w5S55r7uTVQVY4wxpjcKcl0AY4wxA5cFEWOMMb1mQcQYY0yvWRAxxhjTaxZEjDHG9JoFEWOMMb1mQcQYY0yvWRAxxhjTaxZEjDHG9FphrguQaSNHjtTJkyfnuhjGGDOgrFmzplZVq3o6btAHkcmTJ7N69epcF8MYYwYUEXknleMsnWWMMabXLIgYY4zpNQsixhhjem3Qt4kYky3hcJiamhpCoVCui2JMyvx+P+PHj8fr9fbqfAsixvSTmpoaSktLmTx5MiKS6+IY0yNVpa6ujpqaGqZMmdKra1g6y5h+EgqFGDFihAUQM2CICCNGjOhT7dmCiDH9yAKIGWj6+m/WgkgeemVbHZv3Hsp1MYwxpkcWRPLQV3/7Bj95bkuui2EGII/Hw5w5c+KPO+64IyflmDx5MrW1tTl5784efPBBdu/enbX3+9GPfsRDDz0EwDe+8Q1mz57NnDlzmD9/frwcL7zwAuXl5fHP6bbbbkt6rWHDhh2xbcqUKWzatKnDti984QvcddddXZapra2Ns88+m/b29t7eVpcsiOSh5tZ2Wlr7/8M2g18gEGDt2rXxx5IlS3JdpJzrLohEIpF+fa/29nYeeOABrrnmGgC+/OUvs27dOtauXcsHPvCBDsHirLPOin9Ot9xyS8rvcfXVV7Ns2bL462g0yvLly7nqqqu6PMfn83H++efz6KOP9uKuupeTICIiC0Rkk4hsEZEj/pWLyJdEZL2IrBORZ0VkUsK+iIisdR8rslvy7Ai2RQi19+8/bjN0HTx4kOnTp8d/vS5atIif//znAHz605+murqaWbNm8c1vfjN+zuTJk/na177G6aefTnV1Na+++ioXXnghRx99ND/96U8B59f02WefzRVXXMHMmTO54YYbiEajR7z///7v/zJv3jzmzJnDv//7vyf94l6yZAkzZ85k9uzZ3HzzzQBce+213HDDDZx11lkce+yx/P73vwecL/4vf/nLnHLKKcyePZuf/exn8evcddddnHDCCZx44oksWbKE5cuXs3r1aj760Y8yZ84cgsEgkydP5rbbbuPMM8/k8ccf59xzz41PjVRbW0tsrr0HH3yQhQsXcumllzJlyhSWLl3KD37wA+bOnctpp51GfX39Effx3HPPcdJJJ1FY6HR8LSsri+9rbm7udftDbW0tp59+On/4wx9YtGhRhyDyt7/9jcmTJzNp0iQefPBBLr/8chYsWMD06dO59dZb48ctXLiQhx9+uFfv352sd/EVEQ9wD3ABUAOsEpEVqro+4bDXgGpVbRGRTwN3AbEwG1TVOVktdJYFwxGCbRZEBrJb/+8t1u9u7NdrzhxbxjcvndXtMcFgkDlzDv/v8dWvfpWrrrqKpUuXcu211/L5z3+ehoYGPvWpTwFw++23U1lZSSQS4fzzz2fdunXMnj0bgAkTJvDyyy/zxS9+kWuvvZaXXnqJUCjErFmzuOGGGwBYuXIl69evZ9KkSSxYsIAnnniCK6+8Mv7+GzZs4NFHH+Wll17C6/Vy44038vDDD/Ov//qv8WPq6+v57W9/y8aNGxERDhw4EN+3Y8cO/vrXv7J161bOO+88tmzZwkMPPUR5eTmrVq2itbWVM844g/nz57Nx40aefPJJXnnlFYqLi6mvr6eyspKlS5dy9913U11dHb+u3+/nxRdfBIgHxWTefPNNXnvtNUKhEMcccwx33nknr732Gl/84hd56KGH+MIXvtDh+JdeeomTTz65w7avf/3r8TI///zz8e0vv/wyJ554ImPHjuXuu+9m1qzkn+3evXu57LLL+M53vsMFF1wAQEFBAa+//jonnngiy5YtY9GiRfHjV65cyZtvvklxcTGnnHIKl1xyCdXV1Rx//PGsWrWqy3vtrVzUROYBW1R1m6q2AcuAyxMPUNXnVbXFfflPYHyWy5gz4UiUcEQJhY/8RWdMTzqns2IpjgsuuIATTjiBm266ifvvvz9+/GOPPcZJJ53E3Llzeeutt1i//vBvucsuuwyAE044gVNPPZXS0lKqqqrw+/3xL/p58+YxdepUPB4PixYtin8xxzz77LOsWbOGU045hTlz5vDss8+ybdu2DseUlZXh9/v5t3/7N5544gmKi4vj+z7ykY9QUFDAtGnTmDp1Khs3buTPf/4zDz30EHPmzOHUU0+lrq6OzZs385e//IVPfOIT8fMrKyu7/O/UXeon0XnnnRe/7/Lyci699NL4f5MdO3YccfyePXuoquo48e3tt9/Orl27+OhHP8rSpUsBOOmkk3jnnXd4/fXX+exnP8vChQuTvn84HOb888/nrrvuigcQIF4baW9v53e/+x0f/vCH4/suuOACRowYQSAQ4IMf/GD8M/F4PPh8Pg4d6t9OO7kYbDgO2JXwugY4tZvjPwk8lfDaLyKrgXbgDlV9svMJInI9cD3AxIkT+1zgbAqFnRqIpbMGtp5qDNkWjUbZsGEDgUCA+vp6xo8fz/bt27n77rtZtWoVFRUVXHvttR3GCxQVFQHOr97Y89jrWANt5/RM59eqyuLFi/ne977XZdkKCwtZuXIlzz77LMuWLWPp0qU899xzXV5fVfmv//ovLrzwwg77/vSnP6WcLiopKenw/rE0XOfxEp3vO/G/SbJG6kAg0OWYi2uuuYZLLrmEW2+9tUOa6+KLL+bGG2+ktraWkSNHdjinsLCQk08+maeffppzzjknvn3RokXMnz+fc845h9mzZzNq1Kj4vu4+k9bWVvx+f9Ly9VYuaiLJPmVNeqDIx4Bq4PsJmyeqajVwDfAjETn6iIup3qeq1apa3flXQb4LxoKIpbNMP/rhD3/IjBkzeOSRR7juuusIh8M0NjZSUlJCeXk5e/fu5amnnur5Qp2sXLmS7du3E41GefTRRznzzDM77D///PNZvnw5+/btA5zU1TvvdJxhvKmpiYMHD3LxxRfzox/9iLVr18b3Pf7440SjUbZu3cq2bduYPn06F154Iffeey/hcBiAt99+m+bmZubPn88DDzxAS0tL/L0ASktLu/31PXnyZNasWQPA8uXL0/5vkGjGjBls2XK4Z+XmzZvjz1esWMFxxx0HwHvvvYeq87W3cuVKotEoI0aMOOJ6IsIDDzzAxo0bO/S0O/rooxkxYgRLlizpkMoCeOaZZ6ivrycYDPLkk09yxhlnAFBXV0dVVVWvpzfpSi5qIjXAhITX44Ejuk6IyPuBrwPnqGprbLuq7nb/bhORF4C5wNZMFjibQm3uL6J2S2eZ9HVuE1mwYAHXXXcd999/PytXrqS0tJSzzz6b73znO9x6663MnTuXWbNmMXXq1PiXTTpOP/10lixZwhtvvBFvZE80c+ZMvvOd7zB//nyi0Sher5d77rmHSZPifWU4dOgQl19+OaFQCFXlhz/8YXzf9OnTOeecc9i7dy8//elP42mvHTt2cNJJJ6GqVFVV8eSTT7JgwQLWrl1LdXU1Pp+Piy++mO9+97vxBvpAIMDLL798xD3cfPPNfOQjH+FXv/oV73vf+9L+b5Dooosu4uMf/3j89ZIlS9i0aRMFBQVMmjQp3v6yfPly7r33XgoLCwkEAixbtqzLWpTH42HZsmVceumllJWVceONNwJObeSrX/3qEf/NzzzzTD7+8Y+zZcsWrrnmmnhb0PPPP8/FF1/cp/tLSlWz+sAJXNuAKYAPeB2Y1emYWGCY1ml7BVDkPh8JbAZmdvd+J598sg4kG/c06qSv/F5nfOOpXBfFpGn9+vW5LkJWPf/883rJJZdk7PqLFy/Wxx9/PGPXz5SFCxfq22+/nZP3/sUvfqE33XRT0n1XXHGFbty4Mem+ZP92gdWawnd61tNZqtoOfAZ4GtgAPKaqb4nIbSJymXvY94FhwOOduvLOAFaLyOvA8zhtIusZRGLprGA4Eq/uGmMGjjvuuIM9e/bkuhgdtLW1sXDhQqZPn97v15bB/kVVXV2tA2l53Je31rHo5/8EYNN3FlBU6MlxiUyqNmzYwIwZM3JdDGPSluzfroisUaf9uVs2Yj3PxHpnweH2ETNwDPYfZWbw6eu/WQsieSaYGESsm++A4vf7qaurs0BiBgx11xPpS7dfW5QqzySOVE+slZj8N378eGpqati/f3+ui2JMymIrG/aWBZE8k1gTCVoQGVC8Xm+vV4czZqCydFae6dAmYlOfGGPynAWRPJOYzrJJGI0x+c6CSJ6xhnVjzEBiQSTPJAaRVmsTMcbkOQsieSZkDevGmAHEgkieCbZFKCp0PhZrWDfG5DsLInkmGI5QUewDbJyIMSb/WRDJM8FwlOHFXve5BRFjTH6zIJJnQm2ReBCxdJYxJt9ZEMkzwXCEYl8hfm+BpbOMMXnPgkieaWlrJ+D14Pd6LIgYY/KeBZE8EwpH8Xs9BCyIGGMGAAsieSYYjhDwFeD3egj2ok2kPRKlzdZnN8ZkiQWRPBNsi/QpnXXnnzZy1X0vZ6BkxhhzJAsieURVnZqI19PrhvVXdx5gV31LBkpnjDFHsvVE8kirm4YK+ArxF/auJrK9ttnm3DLGZI3VRPJIbOr3gLeAgM+T9jiRAy1t1De30RKO2BKtxpissCCSR2Ij1AM+J52V7oj17bXNAKgertUYY0wmWRDJI7Gg4fd6epXOigURsAWtjDHZYUEkjxxOZ3nw9yKdtSMhiLRYu4gxJgssiOSRUGI6qxc1kW1WEzHGZJkFkTwSbxPxegj40u/iu722GU+BONeyIGKMyQILInkk9sUfaxNpjyrhSGopLVVle20zR1eVAM4cXMYYk2kWRPJIx95ZHiD1han2HWqlpS3CrLHlHa5ljDGZZEEkj4TCHRvWnW2p1URiPbNmjikDLJ1ljMmOnAQREVkgIptEZIuILEmy/0sisl5E1onIsyIyKWHfYhHZ7D4WZ7fkmdWhd1Z8nfXUgkEsiMwa6wSRFgsixpgsyHoQEREPcA9wETATWCQiMzsd9hpQraqzgeXAXe65lcA3gVOBecA3RaQiW2XPtNisvb1JZ22vbcZXWMDRo4a517IgYozJvFzUROYBW1R1m6q2AcuAyxMPUNXnVTU2i+A/gfHu8wuBZ1S1XlUbgGeABVkqd8bFvviLCgsIeNNLZ23b38zkEcUUu2kwS2cZY7IhpQkY3RpAT6KqeiCF48YBuxJe1+DULLrySeCpbs4dl8J7DgghdwZfEYnXRFKtUWyvbeKYUcPiwcfSWcaYbEh1Ft/d7kO6OcYDTEzhWsmukXS2QBH5GFANnJPOuSJyPXA9wMSJqRQpPwTbIgTcmoTfm3qbSHskys76Fi6YeRSFngJ8nvTn3TLGmN5INYhsUNW53R0gIq+leK0aYELC6/E4Aarz9d4PfB04R1VbE849t9O5L3Q+V1XvA+4DqK6uHjDT2cbWEgHSahPZfSBEOKJMHemMEQn4PARtnIgxJgtSbRM5vZ+OAVgFTBORKSLiA64GViQeICJzgZ8Bl6nqvoRdTwPzRaTCbVCf724bFILhSLwGkk46a1ttEwBT3IGGxT6PpbOMMVmRUk1EVUP9cYx7XLuIfAbny98DPKCqb4nIbcBqVV0BfB8YBjwuIgA7VfUyVa0XkW/jBCKA21S1PpX3HQgS01mxv60pNKzHuvdOidVEvB5LZxljsiLtlQ1F5Cuqemdf3lRV/wj8sdO2WxKev7+bcx8AHujL++er2PrqQHycSCrBYHttM6VFhYwo8QGxdJYFEWNM5vUYRETkscSXwBygT0HEJBcMRyj1Ox9JOm0i22ubmVJVgltro9hnNRFjTHakUhNpVNV/i70QkXszWJ4hLRSOMKq0CEgMIqmls06edHjMpd/r4VDIGtaNMZmXSsP67Z1efz0TBTFu7yy3LcRTICl11Q2FI7x7IBhvDwG3JmLpLGNMFvQYRFR1O4CIjHRfD5qG7HyT2CYCUOTteU2RnfUtqNIpiBRaOssYkxXpTHsyKBuz84nTxfdwEAl4PbS2dx8Mtu13emZNHTksvs3vtS6+xpjsSCeIdDda3fSDUEI6C5xg0FNaKta9d/LI4vi2YhtsaIzJknSCyIAZ+T0QhSNRwhHtkM4KeD09Nqxvr21i5LAiSv3eDucFwxFU7SMzxmSW1UTyROKCVDF+b88N6ztqW+LTncQEfB6iCq3tqc0AbIwxvZVSEBGRpcAdGS7LkJa4NG5MkdfTY8P6ttrmDo3qgE0Hb4zJmlRrIpuBu0Vkh4jcKSJzMlmooSjU5i5I1Tmd1U1tojEUprapNT5nVuJ5YAtTGWMyL6Ugoqo/VtXTcaZkrwd+ISIbROQWETk2oyUcIpLVRPzeAkLd1CZ2dJozKyZ2DeuhZYzJtLRWNlTVd1T1Tnda+GuAK4ANGSnZEBNM2ibiIdRNF9/OEy/GFPuciQgsnWWMybS0goiIeEXkUhF5GGe1wbeBD2WkZENM7Au/8ziR7tpEttc2IwITK4s7bLd0ljEmW1JdHvcCYBFwCbASZ13061W1OYNlG1JCSdNZ3Y8T2V7bzLjhgQ6BJ/EaLTZWxBiTYalOBf814NfAzTbtSWZ0nc7qumF9X2MrR5X5j9ge652VygzAxhjTF6kuSnVepgsy1MVqHJ3HibS1R4lEFU/BkcN0DgbDjB1+ZBCJXcMa1o0xmZZum4iIyMdE5Bb39UQRmZeZog0tsZqI33f4I4mlqbqaP+tgMExZwHvE9mLrnWWMyZK0ggjw3zhrqS9yXx8C7unXEg1RyUasB3pYU6QxGKY8SRAJWDrLGJMl6S6Pe6qqniQirwGoaoOI+DJQriEnWe8sv7frJXLbI1EOtbYnDyKWzjLGZEm6NZGwiHhwJ2MUkSrAJmjqB8FwBK9H8HqOTGclq1E0uisXJgsihZ4CfJ4CCyLGmIxLN4j8BPgtMEpEbgdeBL7b76UaglraIkd01e0uiBwMhoHkQcQ5t+cFrYwxpq/STWdtBP4DOB9nVt+FwDH9XaihKBTuuKoh9C2IFPsKbZyIMSbj0q2J/Bzwquo9qroUmAP8Z/8Xa+gJdlqQCrpvWO85iNjqhsaYzEu3JnIlsFxEPgqcCfwrML/fSzUEdV5fHRIa1pMEg1gQGV7cVTqr52nkjTGmr9IKIqq6TUSuBp4EdgHzVTWYkZINMZ3XV4eEdFaScSKxIJJsnAhYTcQYkx2pzp31Bh2Xx60EPMArIoKqzs5E4YaSZG0i3aazWtqArtNZAZ+HQyFrEzHGZFaqNZEPZLQUhmA4wqjSjgGhqJtxIgeDYfzeAooKPUfsAycA7T/U2v8FNcaYBKnOnfVOpgsy1CVvE3GnPekiiHRVCwFLZxljsiPVNdZf7Y9jTNdC4eiRU7r30MW3uyAS8BVaEDHGZFyq6awZIrKum/0ClPdDeYYsp4tvx5ju9RTgKZAu01ndBhHrnWWMyYJUg8hxKRyT8jeWiCwAfozTOH+/qt7Raf/ZwI+A2cDVqro8YV8EeMN9uVNVL0v1ffNZsnQWxIJBsnEi7YxLMg18jJPOakdVETlyGnljjOkPWW8Tcefeuge4AKgBVonIClVdn3DYTuBa4OYklwiq6pz+Kk8+UFWnJpIkiPi9BUlrIo3BMDPHlHV5zYDPQ1ShLRLtsvHdGGP6Kt0R6/1hHrBFVbepahvOUruXJx6gqjtUdR1DZHLHVnf1Qr/vyC/7osLkaalU0lmQfKCiMcb0l1wEkXE4AxVjatxtqfKLyGoR+aeILEx2gIhc7x6zev/+/X0pa1YkW9UwJuDz0NopnRWORGnqYhr4GFuYyhiTDemubPgZEano43smS9Brkm1dmaiq1cA1wI9E5OgjLqZ6n6pWq2p1VVVVb8uZNbF0VXGSmkiydFZjfN6srrORsXm4kqXCjDGmv6RbEzkKpw3jMRFZIL1rsa0BJiS8Hg/sTvVkVd3t/t0GvADM7UUZ8kp8adxkbSJJ0lnxyRe7mDcLLJ1ljMmOtIKIqv4nMA34H5yG780i8t1ktYFurAKmicgUd1XEq4EVqZwoIhUiUuQ+HwmcAazv/qz811M6q3NtoqcZfMGZCh4snWWMyay020RUVYH33Ec7UIEzs+9dKZ7fDnwGeBrYADymqm+JyG0ichmAiJwiIjXAh4Gfichb7ukzgNUi8jrwPHBHp15dA1J8ffUuG9Y7tomkEkRiY04snWWMyaS0ZvEVkc8Bi4Fa4H7gy6oaFpECYDPOglU9UtU/An/stO2WhOercNJcnc/7B3BCOmUeCGJf9F03rKdfEwl4nY82aAtTGWMyKN31REYCH+w8bkRVoyJikzT2UiydlbxNpLuGdV+X17TeWcaYbEg3nVXUOYCIyJ0Aqrqh30o1xAS7SWclW1zqQEsq6SzrnWWMybx0g8gFSbZd1B8FGcpCPaSzkrWJBLwefIVdf3zxIGI1EWNMBqW6KNWngRuBqZ0mYiwFXspEwYaS7npnxdJZiXNg9TRaPfFals4yxmRSqm0ivwaeAr4HLEnYfkhV6/u9VENMS3e9s2JrirQfnio+lSDi9RTg9SSfAdgYY/pLqhMwHgQOAosyW5yhKeTWFoqSpKcS1xRJJ4jEzrV0ljEmk1JdlOpF9+8hEWl0H4dirzNbxMEvNoNvsgkA/EnWWT8YDFOWShDxWRAxxmRWqjWRM92/pZktztDkLEjVxVrp7qDBxB5ajcEw5WN7DiLFvsJ4qswYYzIh3QkYPywipe7z/xSRJ0RkwM9dlWvBtmjSRnVw5s6Cjl11DwTDDO9m3qwYJ51lgw2NMZmTbhffb6jqIRE5E7gQ+CXw0/4v1tDitHck/yj8ndZZD0eitLRFUmsTSTLvljHG9Kd0g0jsG+kS4F5V/R3Q9bBpk5Lu0lmd20RSmfIkxlki14KIMSZz0g0i74rIz4CPAH90Z9TNxcJWg0pX66sD8RpKrCaSThCx3lnGmExLNwB8BGf23QWqegCoBL7c76UaYoIJ3Xc765zOSiuIWDrLGJNhaU3AqKotwBMJr/cAe/q7UENNKBxhVGlR0n3xxaU6BZFUuvhaOssYk2npTgVfBHwImJx4rqre1r/FGlrSaRNpTCudVWjpLGNMRqU7FfzvcEaurwFa+784Q1PG2kR8R867ZYwx/SndIDJeVRdkpCRDWCptIrF0VirTwMcU+wqJRJW2SJSiwuTXN8aYvki3Yf0fIjLoVhbMtVA36ayiwgJEiK9ueDAYptjX/TTwMfEAZCktY0yGpFsTORP4hIhsw0lnCc6y67P7vWRDRDgSJRzRLtNZIoK/0EOo/fA4kVRqIXB4dcNgOMLw/imuMcZ0kG4QsQWo+ll3C1LF+L0F8dpEb4KI9dAyxmRKuumsncBZwGJ3mVwFRvd7qYaQWFuHv4t0FnRcIjfVGXxj54Gls4wxmZNuEPlv4HQOrytyCLinX0s0xITanDRVdzWRgPfwoMHGXqazjDEmE9INIqeq6k1ACEBVG7C5s/ok9gVf3E1NpMjr6TB3lqWzjDH5It0gEhYRD04aCxGpAqLdn2K6E0yxTaS1/XAX31SDiKWzjDGZlm4Q+QnwW2C0iNwOvAh8t99LNYTEvuC7GicChydSbGuPOj2tUq6JOP0mgmFbU8QYkxnpzp31sIisAc53Ny1U1Q39X6yhI947q4eG9cZQ+PBo9RQWpAJLZxljMi+lICIiX+pi10UicpGq/qAfyzSkxL7ge2pYD4WjaU15ApbOMsZkXqo1kdja6tOBU4AV7utLgb/1d6GGklTaRIrccSLpzOALCb2zLIgYYzIkpSCiqrcCiMifgZNU9ZD7+lvA4xkr3RBweJxI181Tfq+H1vZIWjP4Ang9BXg9Qot18TXGZEi6DesTgbaE120408KnRUQWiMgmEdkiIkuS7D9bRF4VkXYRubLTvsUistl9LE73vfNNKMV0VmJNJNUgAk4AspqIMSZT0p325FfAShH5LU433yuAX6ZzAbeL8D3ABUANsEpEVqjq+oTDdgLXAjd3OrcS+CZQ7b7/GvfchjTvI2/EayI9dPENtaffJgJOSsuCiDEmU9LtnXW7iDyFM/UJwCdU9bU033MesEVVtwGIyDLgciAeRFR1h7uv8xiUC4FnVLXe3f8MsAB4JM0y5I1gOILXI3g93aSzCj1Eokptk7OESzpBJOD1WDrLGJMx6dZEUNVXgVf78J7jgF0Jr2uAU/tw7rg+lCXngm1dryUSE+v+u7cxRInP023AOfJcW93QGJM56baJ9IdkS+xpf54rIteLyGoRWb1///60CpdtoXDXqxrGFHljQaQ1rVoIuOksG2xojMmQXASRGmBCwuvxwO7+PFdV71PValWtrqqq6nVBs6G79dVjAt7DNZFUu/cmnms1EWNMpuQiiKwCponIFBHxAVdzeNxJT54G5otIhYhUAPPdbQNWd+urx8TWWd/bGEq7JhLweWzEujEmY9JqExGRIuBDON164+eq6m2pXkNV20XkMzhf/h7gAVV9S0RuA1ar6goROQVnjq4K4FIRuVVVZ6lqvYh8GycQAdwWa2QfqLpbXz3G766P3pDG5IsxTjrLgogxJjPSbVj/HXAQWIOzPG6vqOofgT922nZLwvNVOKmqZOc+ADzQ2/fON6m0iSSmu9KuiVg6yxiTQekGkfGquiAjJRmiguEIo0q7DwyxdBb0IojYOBFjTAal2ybyDxE5ISMlGaJSaRMpKux9TaTY54wTUU21A5wxxqQu3ZrImcC1IrIdJ50lgKrq7H4v2RARCkdTHicCMDzFaeDj53qdgYrhiOIrTNZD2hhjei/dIHJRRkoxhDW1tne7NC50nBIl7S6+sYWp2iL4CnPRGc8YM5il9a2iqu8Aw3GmgL8UGO5uM73Q1Np6fNbPAAAbjklEQVTOwWCYMcP93R6XmO7qTcM6QIsNODTGZEBaQUREPg88DIxyH/8rIp/NRMGGgl31LQBMrCzu9ri+NKzbmiLGmExKN531SeBUVW0GEJE7gZeB/+rvgg0FsSAyoaKHINKHhvWALZFrjMmgdJPkAiR+G0VIPp+VScGuhiAAE3qoiRQUSLw9o7fpLBtwaIzJhHRrIr8AXnHXEwFYyCAa+Jdtu+pbGFZUSEUKPa78hQW0tUfTbli3dJYxJpPSXU/kByLyAk5XX6F364kY1676FsZXBBDpuTLnd7vqpjMNPFg6yxiTWenOnXWnqn6FhPVEEraZNO2sb2HKyJKUjg2kuY5I/Lx4Ost6Zxlj+l+630oXJNlmY0d6QVWpaQj22B4S4y/0pJ3KAiiOjxPpvEikMcb0XUo1ERH5NHAjMFVE1sU2A8OAlzJUtkGttqmNYDjChIpASseXFHl6HNmezOF0ltVEjDH9L9V01q+Bp4DvAUsSth8a6FOx58rO2BiREanVRL7xgZl9SmeFrHeWMSYDUgoiqnoQZwr4RSJyInCWu+vvgAWRXqhpSG2MSMzciRW9eh9fYQGFBWIN68aYjEh3xPrnsBHr/WJnnRNExqcYRPrCVjc0xmRKuuNE/g0bsd4vdjW0UFVa1OP66v0h4PVYOssYkxE2Yj1Hdta3pNyo3lfFVhMxxmRIX0asC3A5NmK9V3bVBzllcu/aOdLl91oQMcZkRl9GrIONWO+VcCTKnoNBJlSOy8r7FfssnWWMyYx0G9aLgONwxocMBy4VkVsyUbDBbPeBIFHteeLF/lLsK7RxIsaYjEg3nfU7nK6+a3CWxzW9sKvenb03Cz2zwEln1TbZx2WM6X/pBpHxqrogIyUZQmIDDSdUZq9h3dJZxphMSLd31j9E5ISMlGQI2dXQQmGBMKY8v3pn7apv4cd/2UwkqlkolTFmMEh17qw3AHWP/4SIbMNJZwmgqjo7c0UcfHbVtzCuIoCnIDu9o0uKCjkQDHMwGO52Uav/fPJN/vr2fs6cNoKTJ1VmpWzGmIEt1ZrIB4BLcWbsPQaY776ObTdp2FXf0uO66v3pirnjCEeiLH1uc5fH/GNLLX99ez8AL2zan62iGWMGuJSCiKq+o6rvAHuBDwE/BH4AfNDdZtKwqyGYlelOYo4fV86HTx7Pg//Ywfba5iP2R6PK957ayLjhAU4cX87zm/ZlrWzGmIEt3TaRh4BZONOcLAVmAL/q70INZk2t7dQ3t2WtUT3m5gunU1To4fY/rD9i3x/e2MMb7x7kSxccy/xZR/Hmu43sOxTKavmMMQNTukFkuqp+UlWfdx/XA8dmomCD1a7YFPBZTGcBjCr1c9N5x/CXDfv4++bD6aq29ijff3oTxx1VysK54zh3ehUAf7WUljEmBekGkddE5LTYCxE5FVuUKi3x7r1ZTGfFXHfmZCZWFvPt36+nPeKsdPjIyp3srG/hKxcdh6dAmDmmjFGlRdYuYoxJSbpB5FScbr47RGQHzgy+54jIGwkrHvZIRBaIyCYR2SIiS5LsLxKRR939r4jIZHf7ZBEJisha9/HTNMufc7viY0SyH0SKCj187eLjeHtvE4+s3MmhUJifPLuZ06eO4NxjnRqIiHDu9Cr+tnl/PNAYY0xX0h1s2OeBhiLiAe7BWa+9BlglIitUNTFZ/0mgQVWPEZGrgTuBq9x9W1V1Tl/LkSs1DUGGFRVSUZz+eun94cJZR3Ha1Ep+8MzbbN3fTF1zG0suOg6Rw92Nz50+isdW1/DqzgPMm2JdfY0xXUurJhLrpdXVI8XLzAO2qOo2VW0DluHMBpzocuCX7vPlwPmS+C03gO2sb2F8RYBc3Y6IcMsHZnEwGObBf+zgktljOHHC8A7HnDltJJ4C4QXrpWWM6UH6i3b33ThgV8LrGndb0mNUtR1nvq4R7r4pIvKaiPxVRM4iCRG5XkRWi8jq/fvzK7ef7TEiycwcW8aieRPxFRbw5fnTj9hf5vdSPamC561dxBjTg1wEkWQ/wTvPs9HVMXuAiao6F/gS8GsRKTviQNX7VLVaVaurqqr6XOD+oqrsamjJSXtIZ9++/Hhe/I/zmDyyJOn+c6ePYsOeRt47aF19jTFdy0UQqQEmJLweD+zu6hgRKQTKgXpVbVXVOgBVXQNsZQB1Md7f1EooHM3aiobdKSgQRpX5u9x/3nFuV9+3LaVljOlaLoLIKmCaiEwRER9wNbCi0zErgMXu8yuB51RVRaTKbZhHRKYC04BtWSp3n8WmgJ84Ivc1kZ5MH13KUWV+6+prjOlWur2z+kxV20XkM8DTgAd4QFXfEpHbgNWqugL4H+BXIrIFqMcJNABnA7eJSDvO+u43qGp9tu+ht3blcIxIukSE846r4vev7yEcieL15OL3hjEm32U9iACo6h+BP3badkvC8xDw4STn/Qb4TcYLmCGxIJLNebP64pxjR/HIyl2seaeB06aO6PkEY8yQYz8vs2hnfQtVpUUEfJ5cFyUlZxwzAq9HbEJGY0yXLIhk0a6GlrxoVE9Vqd/LKZMreWGjtYsYY5KzIJJFu+qDOR8jkq5zp1exae8hXt3ZkOuiGGPykAWRLFm9o57dB4McXTUs10VJy5UnT2BiZTGf+uVqdiRZi8QYM7RZEMmC+uY2PvvIa0ysLGbxGZNzXZy0VJb4+OV181Bg8S9WUtvUmusiGWPyiAWRDItGlf/vsbXUNbVxzzUnUebPzcSLfTFlZAn/s7iavY0hrntwFc2t7bkukjEmT1gQybD7/r6N5zft5xsfmMHx48pzXZxemzuxgnuuOYk33z3ITb9+lbBNE2+MwYJIRq3eUc/3n97EJSeM4WOnTcp1cfrs/Bmjuf2KE3hh036+8pt1bN3fRCgc6fX1VJWDLWE2vtfIC5v28dKWWt589yDvHgjS3NqOaucp1Ywx+SYngw2Hglg7yPiKAN/70Ak5m/q9vy2aN5H3Dob48bObeeLVdwEYVVrE+IoAY4YHKBChPRIlHFHao1HaI0pUFVVQlKg6waOuqY09B0MEuwlCPk8BJUUeAl4Pfp8Hf6GHgM9DsfsoKSqkxFdIcZGHimIfR5X5GV3m56hyP0eV+QfMeBxjBjILIhlwMBjmC4867SBP3PgvA7IdpDtfeP803nfcKLbub6KmIUhNQws1DUHW724EoLBAKPQU4PUIhQVCgQgiIDh/PQUFzBhbxvuOG8VR5X7GlAcYXVZEJKo0tIQ50NLGgWCYhpY2WlojBMMRQu4jGI7Q3NrOvsZWmtvaaW5tp7k1QluS9JrfW+AEIPdRVFjgLAhW4mNEiS/+tzzgZVhRoROUipzgVB7wUjWsiEKb7sWYblkQ6Uet7RF+9fI7LH1+Cwdawtx+xfEDuh2kKyLCiROGH7GYVS41tbbz3sEQextDvHcwxHuNIQ60tBEKR50A1O78PRQKs6u+hbW7DtDQ3EZ7tOuUmadAOKrMz9jhfsYNDzCuIsDkESVMrRrG0VUlDC/2ZfEOjclPFkTSEI0q/7duN8W+QiZUBphQUUxJUWF8+/ef3kRNQ5Czpo3kKwuOG5QBJF8NKyrkmFHDOGZU6uNwVJVDre0cbAnT1OrUaprcms2BYBt7DoTYfSBIzYEgq99p4P/W7SGSEHQqir1MGlESr82UF3spD3id1Fq5G3iGBxhe7B006UxjOrMgkoYXt9Ty+WVrO2wbUeLD7/Xw7oEgM8eU8dB1J3D2sfmzEJbpmohQ5vemnG4MR6Lsqm9he20z2/Y3s622mZ31zbzXGGLje4doDIY5lKT7c8DrYexwP8eMGsb00aVMG13K9KNKmTKyxGZHNgOeBZE0/HNbHZ4C4ZFPncbexhA761uoaWhh/6FWbr7wWC4/cRwFBfaLc7DyegqYWjWMqVXDOH9G8mPCkSgHWsK8dzDEuweC7HYfuxpa2LyviWfW7yVWmSksECpLfFQU+6go8bp/fUwdWcLMsWXMGlNOefHgak8zg48FkTS8sr2e2ePLmTelMtdFMXnK6ymgqrSIqtIiThh/ZDozFI6wbX8zb+89xOZ9h6hraqO+uY2Gljbe3nuIuuY2DrSE48ePGx5g1tgypowsYVxFgLHlTtvM2OEBygMWYEzuWRBJUUtbO+tqDvDJM6fmuihmAPN7PcwcW8bMsWVdHrPvUIgNew6xfncj6/c0sn73QV54ez9t7R17oFUUe+ON/Ee7NaSpVSWMrwhQVGjdm012WBBJ0avvHCAcUU6darUQk1mjSv2MKvVzTkLbWjSq1DW3xVNk7zYE2V7XzNZ9TTy3cT+Pra6JHysCY8sDTB5ZzKQRTlAZVepndFkRo8v8jC71UxYotMZ+0y8siKTole11FAhUT6rIdVHMEFRQIPE02ZwkXasPBsNs3d/Ejtpm3qlr4Z26ZnbUtfDUG3toSEiPxZT4PPGay9SRzt+JlcVUlvgYOWzgLJxmcs+CSIpe2VbP8ePKKR1kAwfN4FAe8HLSxApOmnjkj5xgW4R9h0LsbWxlb6MzlqamIcjW/U2s3tHAitd303mGmYDXQ2WJj9FlRRwzahjTRpUybfQwpo0uZWy532oxJs6CSApC4Qhrdx1g8b8M/PmvzNAT8HmYNKKESSNKku4PhSNsr23m3YYg9c1t1DW3Ud/cSl1zG7sPBHlu474O6TK/t4BRpX5GDvPFa0ejS/1O+8yoEqaMLLE2mSHEgkgKXtt5gLZIlFOnjMh1UYzpd36vhxljypgxpuvG/vrmNrbsa+LtvYfYUdtMbVMr+5ta2V7bzMrt9R1SZgUCEyqLObpqGBMriw8/RhQzoaLYUmWDjAWRFLyyvQ4ROMW69pohqrLEx7wplV12bw+2RdhW28SWfU1s3e80+G/d38Qr2+pobus4yebIYT7GVziBJTbzw5jhAcaUO5NnlhZZo/9AYkEkBa9sq2fGUWXWL9+YLgR8HmaNLWfW2I5jY1SV+uY2djUE2Vnfwi73sbO+hdd2NfCHNzpOJQNQ7PNwVLnfmadsZEm8A8DRVcMYOcxnASbPWBDpQWt7hFd3NvDRU609xJh0iQgjhhUxYljyXmXtkSjvuZNm7jl4+O+eg0G21zbz0pZaWhPGxxT7PEyocGow4yuKmVBZzNSRToAZVxHAYzNGZJ0FkR6sqzlIa3vUxocYkwGFngLGVxQzvqI46f5oVHn3QJBttc1s29/k1mac5Qde3toxVeYrLIgHlEkjit0g4wSbscP91tifIRZEevDKtjoA5k22IGJMthUUCBMqnRrHOZ0mNo2lynbUNXdoi3lz90H+9NZ7HdJkIs5kqRXFPipLfIwY5jwfOzzA0VVOymzSiGILNL1gQaQHr2yv57ijSqkosbUjjMkniamykyd1/JHXHomy91Aru+qdBdN21bew71ArDc3OXGWb3jvkzlnWsVfZxEpnlP/Y4X6OKjvc2D92eIAJlTadTDIWRLoRjkRZvaOBj1SPz3VRjDFpKPQUxNdz6U5Tazvb9zezdX8T2/Y7tZl36pt5a3cjtU2tHY7tPJ3MpMpixg4PMHa4E2RGlfqHZJuMBZFurKs5SDAc4dSpNj7EmMFoWFEhJ4wvTzrjcmt7hH2NrbzXGKKmocWdTqaFHXXN/OnN96hvbutwfGwlzNgqmOOGO7Mtj3dXxBysDf85CSIisgD4MeAB7lfVOzrtLwIeAk4G6oCrVHWHu++rwCeBCPA5VX06U+V8ZbvbHmLjQ4wZcooKPfH2mFOStIkeCoXZ02ndmN0HQrzbEGTl9nreawx1aJfxepz2nSkjSpg8siQ+CHNCZTHjKwL4vQMzVZb1ICIiHuAe4AKgBlglIitUdX3CYZ8EGlT1GBG5GrgTuEpEZgJXA7OAscBfRORYVe04mqmfvLKtnmNGDWPksKJMXN4YM4CV+r2U+r0cO7o06f7Edpl36prZXtvCjtpmdtQ189LWWkLhjlP7jy4r4qgyP1WlfkaVFTGqtIhRpX4qS7wML3Y6Agwv9jK82JtXbTO5qInMA7ao6jYAEVkGXA4kBpHLgW+5z5cDS8UZYXQ5sExVW4HtIrLFvd7L/V3I9kiU1TvqWTh3XH9f2hgzBCS2y5zWKSWuquxvanUHXx4eiLn3UCs1DS28trOBuk7pskQTKgPMGlPO8ePK3EGeZYwcVpSTlVVzEUTGAbsSXtcAp3Z1jKq2i8hBYIS7/Z+dzs3It/xbuxtpbrP2EGNM/xOR+LoxJ3cxjrmtPUpdcyv17mqXDS3O37omZxXMt9yuzImKCgsI+DwEvB78Xg9XzB3H586fltF7yUUQSRYqNcVjUjkXEbkeuB5g4sSJ6ZbPvQa8f8YoTrP2EGNMDvgKCxhTHmBMedc9zBpDYTa4K2AeaAkTCkcIhiME2yKE2qOMLst8Kj4XQaQGmJDwejywu4tjakSkECgH6lM8F1W9D7gPoLq6+oggk4rZ44dz/+JTenOqMcZkRZnfy6lTR+Q0Y1KQg/dcBUwTkSki4sNpKF/R6ZgVwGL3+ZXAc6qq7varRaRIRKYA04CVWSq3McaYTrJeE3HbOD4DPI3TxfcBVX1LRG4DVqvqCuB/gF+5Def1OIEG97jHcBrh24GbMtUzyxhjTM9EO6+LOchUV1fr6tWrc10MY4wZUERkjapW93RcLtJZxhhjBgkLIsYYY3rNgogxxphesyBijDGm1yyIGGOM6bVB3ztLRPYD76Rw6EigNsPFyRa7l/w0WO5lsNwH2L10Z5KqVvV00KAPIqkSkdWpdGcbCOxe8tNguZfBch9g99IfLJ1ljDGm1yyIGGOM6TULIofdl+sC9CO7l/w0WO5lsNwH2L30mbWJGGOM6TWriRhjjOk1CyKAiOwQkTdEZK2IDKjZGkXkARHZJyJvJmyrFJFnRGSz+7cil2VMVRf38i0Redf9bNaKyMW5LGMqRGSCiDwvIhtE5C0R+by7fcB9Lt3cy0D8XPwislJEXnfv5VZ3+xQRecX9XB51l6jIW93cx4Misj3hM5mTlfJYOssJIkC1qg64/uIicjbQBDykqse72+4C6lX1DhFZAlSo6ldyWc5UdHEv3wKaVPXuXJYtHSIyBhijqq+KSCmwBlgIXMsA+1y6uZePMPA+FwFKVLVJRLzAi8DngS8BT6jqMhH5KfC6qt6by7J2p5v7uAH4vaouz2Z5rCYywKnq33DWXEl0OfBL9/kvcf6nz3td3MuAo6p7VPVV9/khYAMwjgH4uXRzLwOOOprcl173ocD7gNgXb95/Lt3cR05YEHEo8GcRWeOuzz7QjVbVPeB8CQCjclyevvqMiKxz0115nwJKJCKTgbnAKwzwz6XTvcAA/FxExCMia4F9wDPAVuCAqra7h9QwAIJk5/tQ1dhncrv7mfxQRDK/wDoWRGLOUNWTgIuAm9y0iskP9wJHA3OAPcD/n9vipE5EhgG/Ab6gqo25Lk9fJLmXAfm5qGpEVecA44F5wIxkh2W3VOnrfB8icjzwVeA44BSgEshKqtSCCKCqu92/+4Df4vzjGsj2urnsWE57X47L02uqutf9HyYK/JwB8tm4uerfAA+r6hPu5gH5uSS7l4H6ucSo6gHgBeA0YLiIxJYKHw/szlW50pVwHwvc1KOqaivwC7L0mQz5ICIiJW6DISJSAswH3uz+rLy3AljsPl8M/C6HZemT2Jeu6woGwGfjNnz+D7BBVX+QsGvAfS5d3csA/VyqRGS4+zwAvB+njed54Er3sLz/XLq4j40JP1AEp10nK5/JkO+dJSJTcWofAIXAr1X19hwWKS0i8ghwLs4MnnuBbwJPAo8BE4GdwIdVNe8brLu4l3NxUiYK7AD+PdaukK9E5Ezg78AbQNTd/DWctoQB9bl0cy+LGHify2ychnMPzg/ox1T1Nvc7YBlOCug14GPur/m81M19PAdUAQKsBW5IaIDPXHmGehAxxhjTe0M+nWWMMab3LIgYY4zpNQsixhhjes2CiDHGmF6zIGKMMabXLIiYQUOc2ZhHpnH8uSLyL5ksUz4QkYg7q+tY93XK3T5F5AURubDTti+IyH+LyNHudTPejdTkLwsiZig7F8jrIJIwkrovgqo6JzYzQ5oeAa7utO1q4BFV3epOvWGGMAsiZkARkckislFEfulONLdcRIoTDvmsiLwqzvowx7nnVIrIk+7x/xSR2e5kgjcAX3R/TZ8lIpNE5Fn3uGdFZKJ7/oMi8hMR+YeIbBORK48omHPcx9x1HtaKyM9ExONubxKR28VZ/+GfIjLa3V4lIr8RkVXu4wx3+7dE5D4R+TPwkIgUi8hjbrkeFWfti2oR+aSI/DDh/T8lIj9IVrYuyjtSRF4WkUvc1192y7FO3DUqcGa3/UBsMj/3v9tYnOnHjbEgYgak6cB9qjobaARuTNhX606meS9ws7vtVuA19/iv4axXsgP4KfBD91f634Gl7r7ZwMPATxKuOwY4E/gAcEfnAonIDOAqnMk85wAR4KPu7hLgn6p6IvA34FPu9h+7738K8CHg/oRLngxcrqrXuPfX4Jbr2+4+cEZZX+bObQXwCZw5k3rkBrI/ALeo6h9EZD4wDWe+pTnAySJytqrWASuBBe6pVwOPqo1SNi4LImYg2qWqL7nP/xfnyz0mNtnhGmCy+/xM4FcAqvocMEJEypNc93Tg1+7zX3W67pOqGlXV9cDoJOeej/PlvkqcKbrPB6a6+9qA3ycp1/uBpe7xK4Cy2DxuwApVDSaUf5lb/jeBde7zZuA5nJrCcYBXVd9IUrbOvMCzwH+o6jPutvnu4zXgVZzZYKe5+xJTWle7r40BnLmijBloOv8KTnwdm/MowuF/35LCNXp6n8S5lJJdT4BfqupXk+wLJ/xyTyxXAXB6QrBwLiQC0NzD+8Xcj1O72kiKtRCgHSeYXQj8NeE9vqeqP0ty/JPAD0TkJCAQW6TKGLCaiBmYJorI6e7zRfScn/8bbmpJRM7FSXk1AoeA0oTj/sHhX9wfTeG6iZ4FrhSRUe77VIrIpB7O+TPwmdgL6XpN7BdxlqNFRGYCJ8R2uIsRTQCuIfUaggLXAceJs0wvwNPAdeKsG4KIjIvdizuJ3wvAA2m8hxkiLIiYgWgDsFhE1uHMvNrTetjfAqrd4+/g8HTs/wdcEWtYBz4HfMI97uM461anxE1z/SfOCpnrcFbNG9P9WXwuVi4RWY/T0J/MfwNV7nW/gpPOOpiw/zHgJVVtSKO8EZyAeZ6I3Kiqf8ZJ5b0sIm/gNKgnBthHgBNx02rGxNgsvmZAcXsH/V5Vj89xUbLG7eXlVdWQiByNU+s5VlXb3P2/x2mgf7aL85tUdVgGy5fR65v8ZjURY/JfMfCiiLyOs/bNp1W1TUSGi8jbOONAkgYQV2PiYMP+EhtsiLP2ixmirCZijDGm16wmYowxptcsiBhjjOk1CyLGGGN6zYKIMcaYXrMgYowxptcsiBhjjOm1/wfCZqhJ/2SGBAAAAABJRU5ErkJggg==\n", - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", @@ -468,7 +447,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -497,7 +476,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -531,7 +510,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -559,7 +538,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -569,7 +548,7 @@ " W_matrix = np.empty((N, M), dtype=object) # to store arbitrary Python objects\n", " for n in range(N): # N and M are \"small\", so looping is OK\n", " for m in range(M):\n", - " W_matrix[n, m] = odl.ScalingOperator(-material_spectra[m][n], space=Y)\n", + " W_matrix[n, m] = odl.ScalingOperator(scalar=-material_spectra[m][n], domain=Y)\n", " \n", " W = odl.ProductSpaceOperator(W_matrix)" ] @@ -583,7 +562,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -603,14 +582,14 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def spectral_projector(ray_trafo, energies, source_spectrum, material_spectra):\n", " ...\n", " \n", - " S_row = [odl.ScalingOperator(energy_weights[i] * source_spectrum[i])\n", + " S_row = [odl.ScalingOperator(scalar=energy_weights[i] * source_spectrum[i], domain=Y)\n", " for i in range(N)]\n", " S = odl.ReductionOperator(*S_row)" ] @@ -624,7 +603,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -645,7 +624,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -674,19 +653,318 @@ " W_matrix = np.empty((N, M), dtype=object) # to store arbitrary Python objects\n", " for n in range(N): # N and M are \"small\", so looping is OK\n", " for m in range(M):\n", - " W_matrix[n, m] = odl.ScalingOperator(-material_spectra[m][n], space=Y)\n", + " W_matrix[n, m] = odl.ScalingOperator(scalar=-material_spectra[m][n], domain=Y)\n", " \n", " W = odl.ProductSpaceOperator(W_matrix)\n", " \n", " exp = odl.ufunc_ops.exp(Y)\n", " U = odl.DiagonalOperator(exp, N)\n", "\n", - " S_row = [odl.ScalingOperator(energy_weights[i] * source_spectrum[i])\n", + " S_row = [odl.ScalingOperator(scalar=energy_weights[i] * source_spectrum[i], domain=Y)\n", " for i in range(N)]\n", " S = odl.ReductionOperator(*S_row) \n", " \n", " return S * U * W * P" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A complete example\n", + "\n", + "Now we add the missing parts like retrieving spectra, defining a ray transform etc.\n", + "\n", + "First some imports, including the [`physdata`](https://github.com/Dih5/physdata) module for retrieving attenuation coefficients from NIST." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! pip install https://github.com/Dih5/physdata/archive/master.zip" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "% matplotlib inline\n", + "import csv\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import physdata.xray\n", + "import odl\n", + "import scipy.interpolate" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We download mass attenuation coefficients (unit $[\\text{cm}^2 / \\text{g}]$) of bone, soft tissue and water from NIST to make a phantom. The density values (unit $[\\text{g} / \\text{cm}^3]$) are from \"the internet\"." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bone_att = np.array(physdata.xray.fetch_coefficients('bone'))[:, :2] # [cm^2/g]\n", + "bone_att[:, 0] *= 1000 # [MV -> kV]\n", + "bone_dens = 1.9 # [g/cm^3]\n", + "\n", + "tissue_att = np.array(physdata.xray.fetch_coefficients('tissue'))[:, :2] # [cm^2/g]\n", + "tissue_att[:, 0] *= 1000 # [MV -> kV]\n", + "tissue_dens = 0.91 # [g/cm^3]\n", + "\n", + "water_att = np.array(physdata.xray.fetch_coefficients('water'))[:, :2] # [cm^2/g]\n", + "water_att[:, 0] *= 1000 # [MV -> kV]\n", + "water_dens = 1.0 # [g/cm^3]\n", + "\n", + "# Alternative: load stored files\n", + "# bone_att = np.load('./data/bone_attenuation.npy')\n", + "# tissue_att = np.load('./data/tissue_attenuation.npy')\n", + "# water_att = np.load('./data/water_attenuation.npy')\n", + "\n", + "plt.figure()\n", + "plt.loglog(bone_att[:, 0], bone_att[:, 1], label='Bone')\n", + "plt.loglog(tissue_att[:, 0], tissue_att[:, 1], label='Tissue')\n", + "plt.loglog(water_att[:, 0], water_att[:, 1], label='Water')\n", + "plt.legend()\n", + "plt.xlabel('photon energy [keV]')\n", + "plt.ylabel('mass attenuation [cm^2/g]')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we read a source spectrum obtained from an [online simulator tool](https://www.oem-xray-components.siemens.com/x-ray-spectra-simulation), using a Tungsten target with 35 kVp and 0.1 mm aluminium filter. The given values are relative photon flux densities (unit $[\\text{mm}^{-2}]$), integrating to 1." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spectrum_file = './data/source_spectrum_0.1mm_AL_1m_air.csv'\n", + "with open(spectrum_file) as f:\n", + " lines = f.readlines()\n", + " lines = [line.replace(',', '.') for line in lines] # Comma decimal marker -> dot\n", + " lines.insert(0, 'keV;photon_density\\n') # Manually add column labels\n", + "\n", + "dialect = csv.Sniffer().sniff(lines[0], delimiters=';')\n", + "reader = csv.DictReader(lines, dialect=dialect)\n", + "\n", + "energies = np.empty(len(lines) - 1)\n", + "photon_dens = np.empty(len(lines) - 1)\n", + "for i, row in enumerate(reader):\n", + " energies[i] = float(row['keV'])\n", + " photon_dens[i] = float(row['photon_density']) # [mm^(-2)]\n", + "\n", + "# Alternative: load stored file\n", + "# source_spectrum = np.load('./data/source_spectrum.npy')\n", + "# voltages = source_spectrum[:, 0]\n", + "# photon_dens = source_spectrum[:, 1]\n", + "\n", + "plt.figure()\n", + "plt.plot(energies, photon_dens, label='source spectrum')\n", + "plt.legend()\n", + "plt.xlabel('photon energy [keV]')\n", + "plt.ylabel('relative photon flux density [1/mm^2]')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we interpolate the material spectra at the given energy values:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x, y = water_att[:, 0], water_att[:, 1]\n", + "water_spectrum = scipy.interpolate.interp1d(x, y)(voltages)\n", + "x, y = bone_att[:, 0], bone_att[:, 1]\n", + "bone_spectrum = scipy.interpolate.interp1d(x, y)(voltages)\n", + "x, y = tissue_att[:, 0], tissue_att[:, 1]\n", + "tissue_spectrum = scipy.interpolate.interp1d(x, y)(voltages)\n", + "\n", + "plt.figure()\n", + "plt.semilogy(voltages, water_spectrum, label='Water')\n", + "plt.semilogy(voltages, bone_spectrum, label='Bone')\n", + "plt.semilogy(voltages, tissue_spectrum, label='Tissue')\n", + "plt.xlabel('photon energy [keV]')\n", + "plt.ylabel('mass attenuation coefficient [cm^2/g]')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we build a phantom from these materials. The phantom will consist of 3 components, that correspond to the separate densities of bone, tissue and water, respectively. We use a simple Shepp-Logan type phantom to keep things short.\n", + "\n", + "ODL provides a tool for creating a phantom from a specification of ellipses. It adds the contribution of all ellipses with provided values, locations, half axes and orientations. We first need to specify in which spatial region and with how many points per axis we want to discretize the phantom. Then we make the separate phantoms my modifying the Shepp-Logan ellipses." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Reconstruction volume: square of 10 cm side length\n", + "reco_space = odl.uniform_discr([-5, -5], [5, 5], (256, 256), dtype='float32')\n", + "\n", + "# Get Shepp-Logan ellipses, a list of lists where each sublist specifies\n", + "# one ellipse\n", + "shepp_logan_ellipses = odl.phantom.transmission.shepp_logan_ellipsoids(\n", + " 2, modified=True)\n", + "\n", + "# Bone: outer shell from ellipses 0 and 1\n", + "bone_indices = [0, 1]\n", + "bone_ellipses = [list(shepp_logan_ellipses[i]) for i in bone_indices]\n", + "bone_ellipses[0][0] = bone_dens\n", + "bone_ellipses[1][0] = -bone_dens\n", + "bone_phantom = odl.phantom.ellipsoid_phantom(reco_space, bone_ellipses)\n", + "\n", + "# Soft tissue: inner part, removing the 2 big tilted ellipses and the\n", + "# 3 small ellipses on the bottom\n", + "tissue_indices = [1, 2, 3, 7, 8, 9]\n", + "tissue_ellipses = [list(shepp_logan_ellipses[i]) for i in tissue_indices]\n", + "tissue_ellipses[0][0] = tissue_dens\n", + "tissue_ellipses[1][0] = -tissue_dens\n", + "tissue_ellipses[2][0] = -tissue_dens\n", + "tissue_ellipses[3][0] = -tissue_dens\n", + "tissue_ellipses[4][0] = -tissue_dens\n", + "tissue_ellipses[5][0] = -tissue_dens\n", + "tissue_phantom = odl.phantom.ellipsoid_phantom(reco_space, tissue_ellipses)\n", + "\n", + "# Water: 3 small ellipses on the bottom\n", + "water_indices = [7, 8, 9]\n", + "water_ellipses = [list(shepp_logan_ellipses[i]) for i in water_indices]\n", + "water_ellipses[0][0] = water_dens\n", + "water_ellipses[1][0] = water_dens\n", + "water_ellipses[2][0] = water_dens\n", + "water_phantom = odl.phantom.ellipsoid_phantom(reco_space, water_ellipses)\n", + "\n", + "bone_phantom.show('Bone phantom', clim=[0, 1.9])\n", + "tissue_phantom.show('Tissue phantom', clim=[0, 1.9])\n", + "water_phantom.show('Water phantom', clim=[0, 1.9])\n", + "_ = (bone_phantom + tissue_phantom + water_phantom).show('Combined phantom', clim=[0, 1.9])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we generate an example geometry (parallel beam 2D) and the corresponding ray transform. We use 180 angles and a detector with 384 pixels (matches the volume pixel size):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "angle_part = odl.uniform_partition(0, np.pi, 180)\n", + "detector_part = odl.uniform_partition(-7.5, 7.5, 384)\n", + "geometry = odl.tomo.Parallel2dGeometry(angle_part, detector_part)\n", + "ray_trafo = odl.tomo.RayTransform(reco_space, geometry)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally we generate some data! To do that we create the spectral projector and call it on our 3-component phantom:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mean_photons_per_pixel = 1e4\n", + "source_spectrum = photon_dens * mean_photons_per_pixel\n", + "material_spectra = [bone_spectrum, tissue_spectrum, water_spectrum]\n", + "fwd_op = spectral_projector(\n", + " ray_trafo, energies, source_spectrum, material_spectra)\n", + "\n", + "# Checking domain and range\n", + "print(repr(fwd_op.domain))\n", + "print(repr(fwd_op.range))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spectral_data = fwd_op([bone_phantom, tissue_phantom, water_phantom])\n", + "# To see something, apply log transform\n", + "log_spectral_data = np.log(mean_photons_per_pixel) - np.log(spectral_data)\n", + "_ = (log_spectral_data).show('Spectral data')\n", + "\n", + "# For comparison the monochromatic data\n", + "log_mono_data = ray_trafo(bone_phantom + tissue_phantom + water_phantom)\n", + "_ = log_mono_data.show('Monochromatic data')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To see the beam hardening effect in action, we compute a simple FBP reconstruction from the data (assuming a linear monochromatic model):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a quick FBP reco operator\n", + "fbp_op = odl.tomo.fbp_op(ray_trafo, filter_type='Hann', frequency_scaling=0.9)\n", + "_ = fbp_op(log_spectral_data).show('FBP from spectral data', clim=[0, 1.9])\n", + "_ = fbp_op(log_mono_data).show('FBP from monochromatic data', clim=[0, 1.9])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Conclusions\n", + "\n", + "- We have a machinery for creating a (presumably) realistic spectral projector from given spectra of source and materials.\n", + "- The difference between spectral and non-spectral data in our example is significant.\n", + "- The implementation is pretty compact and can, in principle, be used directly in numerical solvers since `derivative` and `adjoint` methods are defined automatically.\n", + "\n", + "## Next steps\n", + "\n", + "- The problem is non-linear and quite hard to optimize. A good strategy is required.\n", + "- The implementation works fine but seems to be quite slow. We should identify the bottlenecks!\n", + "- To speed things up, in particular later for the derivative and its adjoint, partial computations can be reused. This requires an implementation of the spectral projector as an `Operator` with manual overrides for `derivative` and `adjoint`.\n", + "- Currently, all energies from 0 to 35 keV are used with equal spacing. We should use less samples in the range after the source intensity peaks.\n", + "- For an energy-discriminating detector we need to handle the case of multiple energy bins, i.e., assign photon energy subsets to detector bins. This changes the operator structure a bit, but not too much." + ] } ], "metadata": {