Explorar o código

interpretation update

spacexerq hai 3 semanas
pai
achega
33f64025f3

+ 3 - 1
services/seq-interp/src/gui/adapters.py

@@ -160,7 +160,9 @@ def validate_timing(hw, seq_data: dict, sync_data: dict) -> List[str]:
     raster = hw.block_duration_raster
     if raster > 0:
         for i, dur in enumerate(sync_data.get("blocks_duration", [])):
-            if dur > 0:
+            # Sub-raster durations are hardware sync/gate ticks, not sequence
+            # blocks - the raster does not apply to them
+            if dur >= raster:
                 rem = dur % raster
                 if 1e-15 < rem < raster - 1e-15:
                     warnings.append(

+ 15 - 5
services/seq-interp/src/interfaces/gradient_exporter.py

@@ -8,11 +8,21 @@ class GradientExporter:
 
     @staticmethod
     def _duplicates_delete(loc_list):
-        new_list = [[0] * 2]
-        for i in range(len(loc_list)):
-            if loc_list[i][0] not in np.transpose(new_list)[0]:
-                new_list.append(loc_list[i])
-        return new_list
+        """
+        Оставляет первую встреченную точку для каждого значения времени
+        (первый столбец), сохраняя исходный порядок; строка с t=0 всегда
+        заменяется ведущей [0, 0].
+        """
+        arr = np.asarray(loc_list, dtype=float)
+        seed = np.zeros((1, 2))
+        if arr.size == 0:
+            return seed
+
+        times = arr[:, 0]
+        _, first_idx = np.unique(times, return_index=True)
+        first_idx.sort()
+        keep = first_idx[times[first_idx] != 0]
+        return np.vstack((seed, arr[keep]))
 
     @staticmethod
     def _gradient_time_convertation(params: dict, time_sample):

+ 9 - 5
services/seq-interp/src/interfaces/pulseq_adapter.py

@@ -61,10 +61,14 @@ class PulseqLoader:
     def _parse_blocks(self, seq) -> list:
         """
         Формирует список блоков из объекта Sequence.
+
+        Работает по сырой таблице событий seq.block_events (столбцы:
+        [?, RF, GX, GY, GZ, ADC, EXT]) — без seq.get_block(), который
+        реконструирует объекты событий и распаковывает shape'ы на каждый
+        блок и на больших последовательностях доминирует по времени.
         """
         blocks = []
-        for block_id in seq.block_events:
-            block = seq.get_block(block_id)
+        for block_id, events in seq.block_events.items():
             duration = seq.block_durations[block_id]
             if duration == 0:
                 self.logger.warning(
@@ -75,12 +79,12 @@ class PulseqLoader:
 
             block_type = []
             has_adc = False
-            if getattr(block, "rf", None) is not None:
+            if events[1]:
                 block_type.append("RF")
-            if getattr(block, "adc", None) is not None:
+            if events[5]:
                 block_type.append("ADC")
                 has_adc = True
-            if any(getattr(block, axis, None) is not None for axis in ("gx", "gy", "gz")):
+            if any(events[2:5]):
                 block_type.append("GRAD")
 
             blocks.append({

+ 18 - 13
services/seq-interp/src/interfaces/rf_exporter.py

@@ -25,11 +25,11 @@ class RFExporter:
         t_rf    = np.asarray(t_rf,    dtype=float)
 
         if len(rf_ampl) == 0:
-            return []
+            return np.empty(0)
 
         rf_ampl_maximum = float(np.max(np.abs(rf_ampl)))
         if rf_ampl_maximum == 0:
-            return []
+            return np.empty(0)
 
         proportional_cf_rf = 127.0 / rf_ampl_maximum
 
@@ -44,12 +44,12 @@ class RFExporter:
         real_dense = np.interp(t_dense, t_rf, rf_ampl.real)
         imag_dense = np.interp(t_dense, t_rf, rf_ampl.imag)
 
-        out_rf_list: list[int] = []
-        for r, i in zip(real_dense, imag_dense):
-            out_rf_list.append(round(r * proportional_cf_rf))
-            out_rf_list.append(round(i * proportional_cf_rf))
-
-        return out_rf_list
+        # Interleaved [re, im, re, im, ...]; np.round is round-half-even,
+        # same as builtin round() used previously.
+        out = np.empty(2 * n_samples)
+        out[0::2] = np.round(real_dense * proportional_cf_rf)
+        out[1::2] = np.round(imag_dense * proportional_cf_rf)
+        return out
 
     def export(self, waveforms: dict, params: dict, output_dir: str):
         rf_raster_local = params["rf_raster_time"]
@@ -63,20 +63,25 @@ class RFExporter:
         else:
             empty_block_time_delay = 0
 
-        rf_out = [0] * int(2 * (empty_block_time_delay // rf_raster_local))
-        rf_out += self._radio_ampl_convertation(
+        head = np.zeros(int(2 * (empty_block_time_delay // rf_raster_local)))
+        body = self._radio_ampl_convertation(
             waveforms["rf"],
             waveforms["t_rf"],
             rf_raster_local,
         )
 
         scale_rf = params.get("scale_rf", 1.0)
-        rf_out = [round(x * scale_rf) for x in rf_out]
+        rf_out = np.round(np.concatenate((head, body)) * scale_rf)
+
+        if len(rf_out) and (rf_out.min() < -128 or rf_out.max() > 127):
+            raise ValueError(
+                f"RF amplitude out of int8 range after scale_rf={scale_rf}: "
+                f"[{rf_out.min():.0f}, {rf_out.max():.0f}]"
+            )
 
         file_path = f"{output_dir}/rf_{rf_raster_local}_raster.bin"
         with open(file_path, "wb") as file_rf:
-            for byte in rf_out:
-                file_rf.write(int(byte).to_bytes(1, byteorder="big", signed=True))
+            file_rf.write(rf_out.astype(np.int8).tobytes())
 
         np.savetxt(f"{output_dir}/rf_time.txt", np.transpose(waveforms["t_rf"]))
         np.savetxt(f"{output_dir}/rf_ampl.txt", np.transpose(waveforms["rf"]))

+ 8 - 1
services/seq-interp/src/interfaces/xml_generator.py

@@ -1,3 +1,5 @@
+from itertools import accumulate
+
 from yattag import Doc, indent
 
 
@@ -24,6 +26,11 @@ class XMLGenerator:
             adc_times_values = []
             adc_times_starts = []
 
+            # Префиксные суммы длительностей: duration_prefix[k] ==
+            # sum(blocks_duration[0:k]). Считаются один раз, иначе суммирование
+            # на каждом ADC-фронте давало O(N^2) на длинных последовательностях.
+            duration_prefix = list(accumulate(blocks_duration, initial=0))
+
             # Массивы синхронизатора имеют ведущий seed-элемент на индексе 0
             # (стартовый блок), значимые блоки идут по индексам 1..nb. Seed
             # представлен заголовочным тегом (*1), поэтому в цикле эмитим
@@ -56,7 +63,7 @@ class XMLGenerator:
                         gate_adc[idx] == 1 and gate_adc[idx - 1] == 0
                     )
                     if is_rising_edge:
-                        adc_times_starts.append(sum(blocks_duration[0:idx]))
+                        adc_times_starts.append(duration_prefix[idx])
                         run_duration = 0
                         run_iter = idx
                         while (run_iter < len(gate_adc)

+ 8 - 8
services/seq-interp/tests/test_adapters.py

@@ -37,14 +37,14 @@ def _make_seq_data(blocks: list, rf=None, t_rf=None,
     return {
         "blocks": blocks,
         "params": {"scale_rf": 1.0},
-        "rf":  rf   or np.array([]),
-        "t_rf": t_rf or np.array([]),
-        "gx":  gx   or np.array([]),
-        "t_gx": t_gx or np.array([]),
-        "gy":  gy   or np.array([]),
-        "t_gy": t_gy or np.array([]),
-        "gz":  gz   or np.array([]),
-        "t_gz": t_gz or np.array([]),
+        "rf":  rf   if rf   is not None else np.array([]),
+        "t_rf": t_rf if t_rf is not None else np.array([]),
+        "gx":  gx   if gx   is not None else np.array([]),
+        "t_gx": t_gx if t_gx is not None else np.array([]),
+        "gy":  gy   if gy   is not None else np.array([]),
+        "t_gy": t_gy if t_gy is not None else np.array([]),
+        "gz":  gz   if gz   is not None else np.array([]),
+        "t_gz": t_gz if t_gz is not None else np.array([]),
     }