Immersed boundary methods have been widely used for simulating flows with complex geometries, as quality boundary-conforming grids are usually difficult to generate for complex geometries, especially when motion and/or deformation is involved. A major task in immersed boundary simulations is to inject the immersed boundary information into the background Cartesian grid, such as the inside/outside status of a grid point with regard to the immersed boundary and the accurate subcell position of the immersed boundary for a grid point next to it. Complex geometries in immersed boundary methods can be conveniently represented with triangulated surfaces placed upon underlying Cartesian grids in a Lagrangian manner. Regular, intuitive implementations using triangulations can be error-prone and/or cumbersome in dealing with robustness issues. In addition, they can be prohibitively expensive for high resolution simulations with complex moving/deforming boundaries. In this paper, a simple, robust, and fast procedure is developed for setting up complex triangulations in immersed boundary simulations. Central to this setup procedure are a ray casting and closest surface point computation algorithms. Several illustrative examples, including high resolution cases with Cartesian grids of up to 2.1 x 10(9) points and triangulations of up to 1.3 x 10(6) surface elements, are performed to demonstrate the robustness and efficiency of our procedure.