@@ -476,6 +476,39 @@ def phys2cvr(
476476 regr_dir = None
477477
478478 if regr_dir is None :
479+ if lag_map :
480+ lag , _ , _ = io .load_nifti_get_mask (lag_map )
481+ if func .shape [:3 ] != lag .shape :
482+ raise ValueError (
483+ f'{ lag_map } and { fname_func } have different sizes!'
484+ )
485+
486+ # Read lag_step and lag_max from file (or try to)
487+ lag = lag * mask
488+
489+ lag_list = np .unique (lag [mask > 0 ])
490+
491+ if lag_step is None :
492+ lag_step = np .unique (np .round (lag_list [1 :] - lag_list [:- 1 ],3 ))
493+ if lag_step .size > 1 :
494+ raise ValueError (
495+ f'phys2cvr found different delta lags in { lag_map } '
496+ )
497+ else :
498+ lag_step = lag_step [0 ]
499+ LGR .warning (
500+ f'phys2cvr detected a delta lag of { lag_step } seconds'
501+ )
502+ else :
503+ LGR .warning (f'Forcing delta lag to be { lag_step } ' )
504+
505+ step = int (lag_step * freq )
506+
507+ if lag_max is None :
508+ lag_max = np .abs (lag_list ).max ()
509+ LGR .warning (f'phys2cvr detected a max lag of { lag_max } seconds' )
510+ else :
511+ LGR .warning (f'Forcing max lag to be { lag_max } ' )
479512 regr , regr_shifts = create_physio_regressor (
480513 func_avg ,
481514 petco2hrf ,
@@ -601,61 +634,27 @@ def phys2cvr(
601634
602635 # If user specified a lag map, use that one to regress things
603636 if lag_map :
604- lag , _ , _ = io .load_nifti_get_mask (lag_map )
605- if func .shape [:3 ] != lag .shape :
606- raise ValueError (
607- f'{ lag_map } and { fname_func } have different sizes!'
608- )
609-
610- # Read lag_step and lag_max from file (or try to)
611- lag = lag * mask
612-
613- lag_list = np .unique (lag )
614-
615- if lag_step is None :
616- lag_step = np .unique (lag_list [1 :] - lag_list [:- 1 ])
617- if lag_step .size > 1 :
618- raise ValueError (
619- f'phys2cvr found different delta lags in { lag_map } '
620- )
621- else :
622- LGR .warning (
623- f'phys2cvr detected a delta lag of { lag_step } seconds'
624- )
625- else :
626- LGR .warning (f'Forcing delta lag to be { lag_step } ' )
627-
628- step = int (lag_step * freq )
629-
630- if lag_max is None :
631- lag_max = np .abs (lag_list ).max ()
632- LGR .warning (f'phys2cvr detected a max lag of { lag_max } seconds' )
633- else :
634- LGR .warning (f'Forcing max lag to be { lag_max } ' )
635-
636- lag_idx = (lag + lag_max ) * freq / step
637-
638- lag_idx_list = np .unique [lag_idx ]
639637
638+ lag_idx = np .round ((lag + lag_max ) * freq / step ).astype (int )
639+ lag_idx_list = np .unique (lag_idx )
640+
640641 # Prepare empty matrices
641642 beta = np .empty_like (lag , dtype = 'float32' )
642643 tstat = np .empty_like (lag , dtype = 'float32' )
643-
644- for i in lag_idx_list :
645- LGR .info (
646- f'Perform L-GLM for lag { lag_list [i ]} ({ i + 1 } of '
647- f'{ len (lag_idx_list )} '
648- )
644+
645+ for index , i in enumerate (lag_idx_list ):
646+ LGR .info (f'Perform L-GLM number { index + 1 } of { len (lag_idx_list )} ' )
649647 regr = regr_shifts [(i * step ), :, np .newaxis ]
650648
651649 x1D = os .path .join (outdir , 'mat' , f'mat_{ i :04g} .1D' )
650+
652651 (beta [lag_idx == i ], tstat [lag_idx == i ], _ ) = stats .regression (
653652 func [lag_idx == i ],
654653 regr ,
655654 denoise_matrix ,
656655 orthogonalised_matrix ,
657656 extra_matrix ,
658- [lag_idx == i ],
657+ [lag_idx == i ][ 0 ] ,
659658 r2model ,
660659 debug ,
661660 x1D ,
@@ -768,7 +767,7 @@ def _main(argv=None):
768767
769768
770769"""
771- Copyright 2021-2025 , Stefano Moia & phys2cvr contributors .
770+ Copyright 2021, Stefano Moia.
772771
773772Licensed under the Apache License, Version 2.0 (the "License");
774773you may not use this file except in compliance with the License.
0 commit comments