MRI-phase-encoding/phase-encoding-simple.py

#!/usr/bin/env python3
# -*- coding: UTF-8 -*-
from scipy.optimize import differential_evolution
import numpy as np
voxel_num = 5
phase_range = np.pi/2
phase_init = 0.0

total_periods = 50
rf_samples_per_period = 7

voxel_amplitudes = np.random.rand(voxel_num)
test_amplitudes = np.zeros(voxel_num)
test_amplitudes = voxel_amplitudes
voxel_phases = np.linspace(0,phase_range, voxel_num)
if len(voxel_amplitudes) != len(voxel_phases):
    print("ERROR! Size of amplitude and phase arrays do not match!")
    raise

omega = 2*np.pi
period = 2*np.pi/omega
time_steps = np.linspace(0, period*total_periods, rf_samples_per_period*total_periods)

def rf_signal(time):
    signal = 0.0
    for i in range(voxel_num):
        signal += voxel_amplitudes[i]*np.sin(
            omega*time + voxel_phases[i] + phase_init
            )
    return signal

def assumed_signal(time, amplitudes):
    signal = 0.0
    for i in range(voxel_num):
        signal += amplitudes[i]*np.sin(
            omega*time + voxel_phases[i] + phase_init
            )
    return signal

def fitness(amplitudes):
    diff = rf_signal(time_steps) - assumed_signal(time_steps, amplitudes)
    return np.sqrt(np.mean(np.square(diff)))

#print(voxel_phases)
print (voxel_amplitudes)

# import matplotlib.pyplot as plt
# plt.plot(time_steps, rf_signal(time_steps))
# plt.show()
#print(fitness(test_amplitudes))

bounds = []
minmax = (0,1)
for i in range(voxel_num):
    bounds.append(minmax)
result = differential_evolution(fitness, bounds, polish=True)
print(result.x, result.fun)
print("Diff")
print((voxel_amplitudes-result.x))
print("percent")
print((voxel_amplitudes-result.x)*100)