diff --git a/src/pe/raytracing/Raytracer.cpp b/src/pe/raytracing/Raytracer.cpp
index ff3cc1c463577649d3a1e3fba05d23966197cf0a..f2dabf8423af22d462fd8b12e095d49518ef1ef9 100644
--- a/src/pe/raytracing/Raytracer.cpp
+++ b/src/pe/raytracing/Raytracer.cpp
@@ -31,7 +31,7 @@ namespace walberla {
 namespace pe {
 namespace raytracing {
    
-void BodyIntersectionInfo_Comparator_MPI_OP( BodyIntersectionInfo_MPI *in, BodyIntersectionInfo_MPI *inout, int *len, MPI_Datatype *dptr) {
+void BodyIntersectionInfo_Comparator_MPI_OP( BodyIntersectionInfo *in, BodyIntersectionInfo *inout, int *len, MPI_Datatype *dptr) {
    for (int i = 0; i < *len; ++i) {
       if (in->bodySystemID != 0 && inout->bodySystemID != 0) {
          WALBERLA_ASSERT(in->imageX == inout->imageX && in->imageY == inout->imageY, "coordinates of infos do not match: " << in->imageX << "/" << in->imageY << " and " << inout->imageX << "/" << inout->imageY);
@@ -88,7 +88,8 @@ Raytracer::Raytracer(const shared_ptr<BlockStorage> forest, BlockDataID storageI
    imageOutputEnabled_(false),
    localImageOutputEnabled_(false),
    filenameTimestepWidth_(5),
-   bodyToShadingParamsFunction_(bodyToShadingParamsFunction) {
+   bodyToShadingParamsFunction_(bodyToShadingParamsFunction),
+   reductionMethod_(MPI_REDUCE) {
    
    setupView_();
    setupFilenameRankWidth_();
@@ -114,7 +115,8 @@ Raytracer::Raytracer(const shared_ptr<BlockStorage> forest, BlockDataID storageI
                      const Config::BlockHandle& config,
                      std::function<ShadingParameters (const BodyID)> bodyToShadingParamsFunction)
    : forest_(forest), storageID_(storageID), globalBodyStorage_(globalBodyStorage),
-   bodyToShadingParamsFunction_(bodyToShadingParamsFunction) {
+   bodyToShadingParamsFunction_(bodyToShadingParamsFunction),
+   reductionMethod_(MPI_REDUCE) {
    WALBERLA_CHECK(config.isValid(), "No valid config passed to raytracer");
    
    pixelsHorizontal_ = config.getParameter<uint16_t>("image_x");
@@ -148,6 +150,13 @@ Raytracer::Raytracer(const shared_ptr<BlockStorage> forest, BlockDataID storageI
    backgroundColor_ = config.getParameter<Color>("backgroundColor", Vec3(0.1, 0.1, 0.1));
 
    blockAABBIntersectionPadding_ = config.getParameter<real_t>("blockAABBIntersectionPadding", real_t(0.0));
+
+   std::string reductionMethod = config.getParameter<std::string>("reductionMethod", "MPI_REDUCE");
+   if (reductionMethod == "MPI_REDUCE") {
+      reductionMethod_ = MPI_REDUCE;
+   } else if (reductionMethod == "MPI_GATHER") {
+      reductionMethod_ = MPI_GATHER;
+   }
       
    setupView_();
    setupFilenameRankWidth_();
@@ -195,13 +204,13 @@ void Raytracer::setupMPI_() {
       MPI_DOUBLE // for color
    };
    MPI_Aint displacements[nblocks];
-   displacements[0] = offsetof(BodyIntersectionInfo_MPI, imageX);
-   displacements[1] = offsetof(BodyIntersectionInfo_MPI, imageY);
-   displacements[2] = offsetof(BodyIntersectionInfo_MPI, bodySystemID);
-   displacements[3] = offsetof(BodyIntersectionInfo_MPI, t);
-   displacements[4] = offsetof(BodyIntersectionInfo_MPI, r);
-   displacements[5] = offsetof(BodyIntersectionInfo_MPI, g);
-   displacements[6] = offsetof(BodyIntersectionInfo_MPI, b);
+   displacements[0] = offsetof(BodyIntersectionInfo, imageX);
+   displacements[1] = offsetof(BodyIntersectionInfo, imageY);
+   displacements[2] = offsetof(BodyIntersectionInfo, bodySystemID);
+   displacements[3] = offsetof(BodyIntersectionInfo, t);
+   displacements[4] = offsetof(BodyIntersectionInfo, r);
+   displacements[5] = offsetof(BodyIntersectionInfo, g);
+   displacements[6] = offsetof(BodyIntersectionInfo, b);
    
    MPI_Datatype tmp_type;
    MPI_Type_create_struct(nblocks, blocklengths, displacements, types, &tmp_type);
@@ -343,7 +352,7 @@ void Raytracer::writeImageBufferToFile(const std::vector<Color>& imageBuffer, co
    }
 }
    
-void Raytracer::syncImageUsingMPIReduce(std::vector<BodyIntersectionInfo_MPI>& intersectionsBuffer, WcTimingTree* tt) {
+void Raytracer::syncImageUsingMPIReduce(std::vector<BodyIntersectionInfo>& intersectionsBuffer, WcTimingTree* tt) {
    WALBERLA_MPI_BARRIER();
    if (tt != NULL) tt->start("Reduction");
    int rank = mpi::MPIManager::instance()->rank();
@@ -363,6 +372,43 @@ void Raytracer::syncImageUsingMPIReduce(std::vector<BodyIntersectionInfo_MPI>& i
    WALBERLA_MPI_BARRIER();
    if (tt != NULL) tt->stop("Reduction");
 }
+   
+void Raytracer::syncImageUsingMPIGather(std::vector<BodyIntersectionInfo>& intersections, std::vector<BodyIntersectionInfo>& intersectionsBuffer, WcTimingTree* tt) {
+   WALBERLA_MPI_BARRIER();
+   if (tt != NULL) tt->start("Reduction");
+   
+   mpi::SendBuffer sendBuffer;
+   for (auto& info: intersections) {
+      sendBuffer << info.imageX << info.imageY
+      << info.bodySystemID << info.t
+      << info.r << info.g << info.b;
+   }
+
+   mpi::RecvBuffer recvBuffer;
+   mpi::gathervBuffer(sendBuffer, recvBuffer, 0);
+   
+   WALBERLA_ROOT_SECTION() {
+      BodyIntersectionInfo info;
+      while (!recvBuffer.isEmpty()) {
+         recvBuffer >> info.imageX;
+         recvBuffer >> info.imageY;
+         recvBuffer >> info.bodySystemID;
+         recvBuffer >> info.t;
+         recvBuffer >> info.r;
+         recvBuffer >> info.g;
+         recvBuffer >> info.b;
+         
+         size_t i = coordinateToArrayIndex(info.imageX, info.imageY);
+         
+         if (intersectionsBuffer[i].bodySystemID == 0 || info.t < intersectionsBuffer[i].t) {
+            intersectionsBuffer[i] = info;
+         }
+      }
+   }
+   
+   WALBERLA_MPI_BARRIER();
+   if (tt != NULL) tt->stop("Reduction");
+}
 }
 }
 }
diff --git a/src/pe/raytracing/Raytracer.h b/src/pe/raytracing/Raytracer.h
index 853c72cc7e3141a4c4fb44a4d46edb854ca9b99c..2b57c0a0230a0e92c5b33a7ca408f2cded12952e 100644
--- a/src/pe/raytracing/Raytracer.h
+++ b/src/pe/raytracing/Raytracer.h
@@ -49,14 +49,6 @@ namespace raytracing {
 /*!\brief Contains information about a ray-body intersection.
  */
 struct BodyIntersectionInfo {
-   size_t imageX;                //!< Viewing plane x pixel coordinate the ray belongs to.
-   size_t imageY;                //!< Viewing plane y pixel coordinate the ray belongs to.
-   walberla::id_t bodySystemID;  //!< System ID of body which was hit.
-   real_t t;                     //!< Distance from camera to intersection point on body.
-   Vec3 color;                   //!< Color computed for the pixel.
-};
-
-struct BodyIntersectionInfo_MPI {
    unsigned int imageX;          //!< Viewing plane x pixel coordinate the ray belongs to.   -> MPI_UNSIGNED
    unsigned int imageY;          //!< Viewing plane y pixel coordinate the ray belongs to.   -> MPI_UNSIGNED
    walberla::id_t bodySystemID;  //!< System ID of body which was hit.                       -> MPI_UNSIGNED_LONG_LONG
@@ -67,6 +59,8 @@ struct BodyIntersectionInfo_MPI {
 };
    
 class Raytracer {
+   enum ReductionMethod { MPI_REDUCE, MPI_GATHER };
+   
 public:
    /*!\name Constructors */
    //@{
@@ -118,6 +112,7 @@ private:
    std::function<ShadingParameters (const BodyID)> bodyToShadingParamsFunction_; /*!< Function which returns a 
                                                                                   * ShadingParameters struct
                                                                                   * given the specified body. */
+   ReductionMethod reductionMethod_; //!< Reduction method used for assembling the image from all processes.
    //@}
    
    /*!\name Member variables for raytracing geometry */
@@ -164,6 +159,7 @@ public:
    inline void setLocalImageOutputEnabled(const bool enabled);
    inline void setImageOutputDirectory(const std::string& path);
    inline void setFilenameTimestepWidth(int8_t width);
+   inline void setReductionMethod(ReductionMethod reductionMethod);
    //@}
    
    /*!\name Functions */
@@ -182,7 +178,8 @@ private:
    void writeImageBufferToFile(const std::vector<Color>& imageBuffer, size_t timestep, bool isGlobalImage = false) const;
    void writeImageBufferToFile(const std::vector<Color>& imageBuffer, const std::string& fileName) const;
    
-   void syncImageUsingMPIReduce(std::vector<BodyIntersectionInfo_MPI>& intersectionsBuffer, WcTimingTree* tt = NULL);
+   void syncImageUsingMPIReduce(std::vector<BodyIntersectionInfo>& intersectionsBuffer, WcTimingTree* tt = NULL);
+   void syncImageUsingMPIGather(std::vector<BodyIntersectionInfo>& intersections, std::vector<BodyIntersectionInfo>& intersectionsBuffer, WcTimingTree* tt);
    
    inline bool isPlaneVisible(const PlaneID plane, const Ray& ray) const;
    inline size_t coordinateToArrayIndex(size_t x, size_t y) const;
@@ -362,7 +359,11 @@ inline void Raytracer::setImageOutputDirectory(const std::string& path) {
 inline void Raytracer::setFilenameTimestepWidth(int8_t width) {
    filenameTimestepWidth_ = width;
 }
-   
+
+inline void Raytracer::setReductionMethod(ReductionMethod reductionMethod) {
+   reductionMethod_ = reductionMethod;
+}
+
 /*!\brief Checks if a plane should get rendered.
  * \param plane Plane to check for visibility.
  * \param ray Ray which is intersected with plane.
@@ -395,8 +396,6 @@ void Raytracer::rayTrace(const size_t timestep, WcTimingTree* tt) {
    if (tt != NULL) tt->start("Raytracing");
    if (tt != NULL) tt->start("Init");
 
-   bool useMPIReduce = true;
-
    real_t inf = std::numeric_limits<real_t>::max();
 
    int numProcesses = mpi::MPIManager::instance()->numProcesses();
@@ -406,7 +405,7 @@ void Raytracer::rayTrace(const size_t timestep, WcTimingTree* tt) {
    std::vector<Color> imageBuffer(pixelsVertical_ * pixelsHorizontal_);
    std::vector<BodyIntersectionInfo> intersections; // contains for each pixel information about an intersection, if existent
    
-   std::vector<BodyIntersectionInfo_MPI> intersectionsBuffer(pixelsVertical_ * pixelsHorizontal_);
+   std::vector<BodyIntersectionInfo> intersectionsBuffer(pixelsVertical_ * pixelsHorizontal_);
 
    
    real_t t, t_closest;
@@ -482,26 +481,16 @@ void Raytracer::rayTrace(const size_t timestep, WcTimingTree* tt) {
             Color color = getColor(body_closest, ray, t_closest, n_closest);
             imageBuffer[coordinateToArrayIndex(x, y)] = color;
             BodyIntersectionInfo intersectionInfo = {
-               x,
-               y,
+               uint32_t(x),
+               uint32_t(y),
                body_closest->getSystemID(),
                t_closest,
-               color
+               color[0],
+               color[1],
+               color[2]
             };
             
-            if (useMPIReduce) {
-               BodyIntersectionInfo_MPI mpi_intersectionInfo = {
-                  static_cast<unsigned int>(x),
-                  static_cast<unsigned int>(y),
-                  body_closest->getSystemID(),
-                  t_closest,
-                  color[0],
-                  color[1],
-                  color[2]
-               };
-               intersectionsBuffer[coordinateToArrayIndex(x, y)] = mpi_intersectionInfo;
-            }
-            
+            intersectionsBuffer[coordinateToArrayIndex(x, y)] = intersectionInfo;
             intersections.push_back(intersectionInfo);
          }
          
@@ -510,51 +499,15 @@ void Raytracer::rayTrace(const size_t timestep, WcTimingTree* tt) {
    }
    if (tt != NULL) tt->stop("Intersection Testing");
    
-   if (useMPIReduce) {
-      syncImageUsingMPIReduce(intersectionsBuffer, tt);
+   switch(reductionMethod_) {
+      case MPI_REDUCE:
+         syncImageUsingMPIReduce(intersectionsBuffer, tt);
+         break;
+      case MPI_GATHER:
+         syncImageUsingMPIGather(intersections, intersectionsBuffer, tt);
+         break;
    }
    
-   /*mpi::SendBuffer sendBuffer;
-   for (auto& info: intersections) {
-      sendBuffer << info.imageX << info.imageY
-      << info.bodySystemID << info.t
-      << info.color[0] << info.color[1] << info.color[2];
-   }
-   int gatheredIntersectionCount = 0;
-   std::vector<real_t> fullTBuffer(pixelsHorizontal_ * pixelsVertical_, inf);
-   std::vector<Color> fullImageBuffer(pixelsHorizontal_ * pixelsVertical_, backgroundColor_);
-   mpi::RecvBuffer recvBuffer;
-   
-   mpi::gathervBuffer(sendBuffer, recvBuffer, 0);
-   //mpi::allGathervBuffer(sendBuffer, recvBuffer);
-   
-   WALBERLA_ROOT_SECTION() {
-      BodyIntersectionInfo info;
-      while (!recvBuffer.isEmpty()) {
-         recvBuffer >> info.imageX;
-         recvBuffer >> info.imageY;
-         recvBuffer >> info.bodySystemID;
-         recvBuffer >> info.t;
-         recvBuffer >> info.color[0];
-         recvBuffer >> info.color[1];
-         recvBuffer >> info.color[2];
-         
-         size_t i = coordinateToArrayIndex(info.imageX, info.imageY);
-         real_t currentFullTBufferT = fullTBuffer[i];
-         if (currentFullTBufferT > info.t) {
-            fullTBuffer[i] = info.t;
-            fullImageBuffer[i] = info.color;
-         }
-         
-         gatheredIntersectionCount++;
-      }
-   }
-   if (tt != NULL) tt->stop("Reduction");
-
-   //WALBERLA_LOG_INFO("#particles visible: " << visibleBodyIDs.size());
-   WALBERLA_LOG_INFO_ON_ROOT("#gathered intersections: " << gatheredIntersectionCount);
-   */
-   
    if (tt != NULL) tt->start("Output");
    if (getImageOutputEnabled()) {
       if (getLocalImageOutputEnabled()) {
@@ -563,10 +516,8 @@ void Raytracer::rayTrace(const size_t timestep, WcTimingTree* tt) {
       WALBERLA_ROOT_SECTION() {
          std::vector<Color> fullImageBuffer(pixelsHorizontal_ * pixelsVertical_, backgroundColor_);
 
-         if (useMPIReduce) {
-            for (auto& info: intersectionsBuffer) {
-               fullImageBuffer[coordinateToArrayIndex(info.imageX, info.imageY)] = Color(info.r, info.g, info.b);
-            }
+         for (auto& info: intersectionsBuffer) {
+            fullImageBuffer[coordinateToArrayIndex(info.imageX, info.imageY)] = Color(info.r, info.g, info.b);
          }
          
          writeImageBufferToFile(fullImageBuffer, timestep, true);
@@ -578,10 +529,8 @@ void Raytracer::rayTrace(const size_t timestep, WcTimingTree* tt) {
       WALBERLA_ROOT_SECTION() {
          std::vector<real_t> fullTBuffer(pixelsHorizontal_ * pixelsVertical_, inf);
 
-         if (useMPIReduce) {
-            for (auto& info: intersectionsBuffer) {
-               fullTBuffer[coordinateToArrayIndex(info.imageX, info.imageY)] = info.t;
-            }
+         for (auto& info: intersectionsBuffer) {
+            fullTBuffer[coordinateToArrayIndex(info.imageX, info.imageY)] = info.t;
          }
          
          writeTBufferToFile(fullTBuffer, timestep, true);