Index: trunk/psModules/src/objects/Makefile.am
===================================================================
--- trunk/psModules/src/objects/Makefile.am	(revision 26076)
+++ trunk/psModules/src/objects/Makefile.am	(revision 26162)
@@ -60,4 +60,5 @@
 	pmSourceMatch.c \
 	pmDetEff.c \
+	pmSourceGroups.c \
 	models/pmModel_GAUSS.c \
 	models/pmModel_PGAUSS.c \
@@ -96,4 +97,5 @@
 	pmSourceMatch.h \
 	pmDetEff.h \
+	pmSourceGroups.h \
 	models/pmModel_GAUSS.h \
 	models/pmModel_PGAUSS.h \
Index: trunk/psModules/src/objects/pmSourceGroups.c
===================================================================
--- trunk/psModules/src/objects/pmSourceGroups.c	(revision 26162)
+++ trunk/psModules/src/objects/pmSourceGroups.c	(revision 26162)
@@ -0,0 +1,117 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+
+#include "pmFPA.h"
+#include "pmSource.h"
+#include "pmSourceGroups.h"
+
+// the strategy here is to divide the image into 2x2 blocks of cells and cycle through
+// the four discontiguous sets of cells, threading all within a set and blocking between
+// sets
+
+// we divide the image region into 2*2 blocks of size Nx*Ny, the image will have
+// Cx*Cy blocks so that (2Nx)Cx = numCols, (2Ny)Cy = numRows.  We want to choose Cx and
+// Cy so that (2Nx)Cx * (2Ny)Cy = 4 * NFILL * nThreads -- each of the four sets of cells
+// has enough cells to allow NFILL cells for each thread (to better distribute heavy and
+// light load cells
+
+// the array runs from readout->image->col0 to readout->image->col0 + readout->image->numCols
+
+
+static void sourceGroupsFree(pmSourceGroups *groups)
+{
+    psFree(groups->groups);
+    return;
+}
+
+
+bool pmSourceGroupsCoordToCell(int *group, // Group number, returned
+                               int *cell,  // Cell number, returned
+                               float x, float y, // Coordinates
+                               const pmSourceGroups *groups // Groups
+                               )
+{
+    // XXX need to handle edges
+    int ix = (x - groups->Xo) / (2 * groups->Nx);
+    ix = PS_MAX(0, PS_MIN(ix, groups->Cx - 1));
+
+    int iy = (y - groups->Yo) / (2 * groups->Ny);
+    iy = PS_MAX(0, PS_MIN(iy, groups->Cy - 1));
+
+    int jx = (((int)(x - groups->Xo)) % (2 * groups->Nx)) / groups->Nx;
+    jx = PS_MAX(0, PS_MIN(jx, groups->Nx - 1));
+
+    int jy = (((int)(y - groups->Yo)) % (2 * groups->Ny)) / groups->Ny;
+    jy = PS_MAX(0, PS_MIN(jy, groups->Ny - 1));
+
+    *group = jx + 2 * jy;
+    *cell  = ix + groups->Cx * iy;
+
+    return true;
+}
+
+
+pmSourceGroups *pmSourceGroupsAlloc(const pmReadout *readout, const psArray *sources, int nThreads)
+{
+    pmSourceGroups *groups = psAlloc(sizeof(pmSourceGroups)); // Groups, to return
+    psMemSetDeallocator(groups, (psFreeFunc)sourceGroupsFree);
+
+    int nCells = nThreads * 2*2;        // number of cells in a single set
+    int C = sqrt(nCells) + 0.5;
+    int Cx, Cy;                         // Number of cells in x and y
+
+    // we need to assign Cx and Cy based on the dimensionality of the image
+    // crude way to find most evenly balanced factors of nCells:
+    if (C <= 1) {
+        Cx = Cy = 1;
+    } else {
+        for (int i = C; i >= 1; i--) {
+            int C1 = nCells / C;
+            int C2 = nCells / C1;
+            if (C1*C2 != nCells) continue;
+
+            if (readout->image->numRows > readout->image->numCols) {
+                Cx = PS_MAX(C1, C2);
+                Cy = PS_MIN(C1, C2);
+            } else {
+                Cx = PS_MAX(C1, C2);
+                Cy = PS_MIN(C1, C2);
+            }
+        }
+    }
+
+    groups->Cx = Cx;
+    groups->Cy = Cy;
+    groups->Xo = readout->image->col0, groups->Yo = readout->image->row0; // Offsets
+    groups->Nx = readout->image->numCols / (Cx*2);
+    groups->Ny = readout->image->numRows / (Cy*2);
+    groups->groups = psArrayAlloc(4);
+
+    long numSources = sources->n;       // Number of input sources
+
+    // Populate the groups
+    for (int i = 0; i < groups->groups->n; i++) {
+        psArray *cells = psArrayAlloc(Cx * Cy); // Cells within a group
+        groups->groups->data[i] = cells;
+        for (int j = 0; j < cells->n; j++) {
+            cells->data[j] = psArrayAllocEmpty(numSources / (Cx * Cy));
+        }
+    }
+
+    // Populate the cells
+    for (int i = 0; i < numSources; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        int group = 0, cell = 0;        // Group and cell index for source
+        pmSourceGroupsCoordToCell(&group, &cell, source->peak->xf, source->peak->yf, groups);
+
+        psArray *cells = groups->groups->data[group]; // Cells for group
+        psArray *cellSources = cells->data[cell];     // Array of sources for cell within group
+        psArrayAdd(cellSources, cellSources->nalloc, source);
+    }
+
+    return groups;
+}
Index: trunk/psModules/src/objects/pmSourceGroups.h
===================================================================
--- trunk/psModules/src/objects/pmSourceGroups.h	(revision 26162)
+++ trunk/psModules/src/objects/pmSourceGroups.h	(revision 26162)
@@ -0,0 +1,39 @@
+#ifndef PM_SOURCE_GROUPS_H
+#define PM_SOURCE_GROUPS_H
+
+#include <pslib.h>
+
+#include "pmFPA.h"
+
+/// Groups of sources
+///
+/// We divide up the sources for threading.
+typedef struct {
+    int Xo, Yo;                         // Offset
+    int Nx, Ny;                         // Size of cells
+    int Cx, Cy;                         // Number of cells in x and y
+    psArray *groups;                    // Cell groups
+} pmSourceGroups;
+
+/// Allocator
+///
+/// Allocates the source groups
+pmSourceGroups *pmSourceGroupsAlloc(
+    const pmReadout *readout,           // Readout on which the sources are defined
+    const psArray *sources,             // Sources of interest
+    int nThreads                        // Number of threads
+    );
+
+
+/// Return the group and cell indices given x,y coordinates
+bool pmSourceGroupsCoordToCell(
+    int *group,                         // Group number, returned
+    int *cell,                          // Cell number, returned
+    float x, float y,                   // Coordinates
+    const pmSourceGroups *groups        // Groups
+    );
+
+
+
+
+#endif