docs: def _create_features(self, constraints):
      """Creates all the matrices needed to compute the output warped grids."""
      affine_warp_constraints = constraints
      if not isinstance(affine_warp_constraints, AffineWarpConstraints):
        affine_warp_constraints = AffineWarpConstraints(affine_warp_constraints)
      mask = affine_warp_constraints.mask
      psi = _create_affine_features(output_shape=self._output_shape,
                                    source_shape=self._source_shape)
      scales = [(x - 1.0) * .5 for x in reversed(self._source_shape)]
      offsets = scales
      # Transforming a point x's i-th coordinate via an affine transformation
      # is performed via the following dot product:
      #
      #  x_i' = s_i * (T_i * x) + t_i                                          (1)
      #
      # where Ti is the i-th row of an affine matrix, and the scalars s_i and t_i
      # define a decentering and global scaling into the source space.
      # In the AffineGridWarper some of the entries of Ti are provided via the
      # input, some others are instead fixed, according to the constraints
      # assigned in the constructor.
      # In create_features the internal dot product (1) is accordingly broken down
      # into two parts:
      #
      # x_i' = Ti[uncon_i] * x[uncon_i, :] + offset(con_var)                   (2)
      #
      # i.e. the sum of the dot product of the free parameters (coming
      # from the input) indexed by uncond_i and an offset obtained by
      # precomputing the fixed part of (1) according to the constraints.
      # This step is implemented by analyzing row by row the constraints matrix
      # and saving into a list the x[uncon_i] and offset(con_var) data matrices
      # for each output dimension.
      features = []
      for row, scale in zip(mask, scales):
        x_i = np.array([x for x, is_active in zip(psi, row) if is_active])
        features.append(x_i * scale if len(x_i) else None)
  
      for row_i, row in enumerate(mask):
        x_i = None
        s = scales[row_i]
        for i, is_active in enumerate(row):
          if is_active:
            continue
  
          # In principle a whole row of the affine matrix can be fully
          # constrained. In that case the corresponding dot product between input
          # parameters and grid coordinates doesn't need to be implemented in the
          # computation graph since it can be precomputed.
          # When a whole row if constrained, x_i - which is initialized to
          # None - will still be None at the end do the loop when it is appended
          # to the features list; this value is then used to detect this setup
          # in the build function where the graph is assembled.
          if x_i is None:
            x_i = np.array(psi[i]) * affine_warp_constraints[row_i][i] * s
          else:
            x_i += np.array(psi[i]) * affine_warp_constraints[row_i][i] * s
        features.append(x_i)
  
      features += offsets
      return features