@@ -263,8 +263,31 @@ def trace_cubes(self, beam: Union[dict[str, Any], Beam]) -> list[sitk.Image]:
263263        for  i , cube  in  enumerate (rad_depth_cubes ):
264264            rel_depths  =  lengths  *  rho [i ]
265265            rel_depths  =  np .cumsum (rel_depths , axis = 1 ) -  rel_depths  /  2.0 
266-             ix_assign  =  np .unravel_index (ix [ix_remember_from_tracing ], cube .shape , order = "F" )
267-             cube [ix_assign ] =  rel_depths [ix_remember_from_tracing ]
266+ 
267+             try :
268+                 ix_assign  =  np .unravel_index (ix [ix_remember_from_tracing ], cube .shape , order = "F" )
269+             except  (ValueError , IndexError ):
270+                 logger .error (
271+                     "Error in unraveling indices from raytracing. Trying to recover..." ,
272+                     exc_info = True ,
273+                 )
274+                 tmp_ix  =  ix [ix_remember_from_tracing ]
275+                 rel_depths  =  rel_depths [ix_remember_from_tracing ]
276+ 
277+                 wrong_values  =  np .logical_or (tmp_ix  <  0 , tmp_ix  >=  cube .size )
278+                 tmp_ix  =  tmp_ix [~ wrong_values ]
279+                 rel_depths  =  rel_depths [~ wrong_values ]
280+ 
281+                 # Remove the wrong values 
282+                 ix_assign  =  np .unravel_index (tmp_ix , cube .shape , order = "F" )
283+                 logger .info (
284+                     "Recovered %d indices for radiological depth cube" ,
285+                     np .count_nonzero (wrong_values ),
286+                 )
287+                 cube [ix_assign ] =  rel_depths 
288+             else :
289+                 cube [ix_assign ] =  rel_depths [ix_remember_from_tracing ]
290+ 
268291            rad_depth_cubes [i ] =  sitk .GetImageFromArray (cube )
269292            rad_depth_cubes [i ].CopyInformation (self .cubes [i ])
270293
0 commit comments