diff --git a/src/library.cpp b/src/library.cpp
index ca3276ee8ed59e486ce71924bf3ced417c147c26..892818ff9ca873509226b13ab4b67ebd9b4e4b3a 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -773,7 +773,9 @@ void lammps_gather_atoms(void *ptr, char *name,
     if (type == 0) {
       int *vector = NULL;
       int **array = NULL;
-      if (count == 1) vector = (int *) vptr;
+      const int imgunpack = (count == 3) && (strcmp(name,"image") == 0);
+
+      if ((count == 1) || imgunpack) vector = (int *) vptr;
       else array = (int **) vptr;
 
       int *copy;
@@ -786,11 +788,19 @@ void lammps_gather_atoms(void *ptr, char *name,
       if (count == 1)
         for (i = 0; i < nlocal; i++)
           copy[tag[i]-1] = vector[i];
-      else
+      else if (imgunpack) {
+        for (i = 0; i < nlocal; i++) {
+          offset = count*(tag[i]-1);
+          const int image = vector[i];
+          copy[offset++] = (image & IMGMASK) - IMGMAX;
+          copy[offset++] = ((image >> IMGBITS) & IMGMASK) - IMGMAX;
+          copy[offset++] = ((image >> IMG2BITS) & IMGMASK) - IMGMAX;
+        }
+      } else
         for (i = 0; i < nlocal; i++) {
           offset = count*(tag[i]-1);
           for (j = 0; j < count; j++)
-            copy[offset++] = array[i][0];
+            copy[offset++] = array[i][j];
         }
       
       MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);
@@ -874,7 +884,9 @@ void lammps_scatter_atoms(void *ptr, char *name,
     if (type == 0) {
       int *vector = NULL;
       int **array = NULL;
-      if (count == 1) vector = (int *) vptr;
+      const int imgpack = (count == 3) && (strcmp(name,"image") == 0);
+
+      if ((count == 1) || imgpack) vector = (int *) vptr;
       else array = (int **) vptr;
       int *dptr = (int *) data;
 
@@ -882,6 +894,15 @@ void lammps_scatter_atoms(void *ptr, char *name,
         for (i = 0; i < natoms; i++)
           if ((m = lmp->atom->map(i+1)) >= 0)
             vector[m] = dptr[i];
+      } else if (imgpack) {
+        for (i = 0; i < natoms; i++)
+          if ((m = lmp->atom->map(i+1)) >= 0) {
+            offset = count*i;
+            int image = dptr[offset++] + IMGMAX;
+            image += (dptr[offset++] + IMGMAX) << IMGBITS;
+            image += (dptr[offset++] + IMGMAX) << IMG2BITS;
+            vector[m] = image;
+          }
       } else {
         for (i = 0; i < natoms; i++)
           if ((m = lmp->atom->map(i+1)) >= 0) {