diff --git a/src/cristallina/analysis.py b/src/cristallina/analysis.py
index 3a773b8..36151ef 100644
--- a/src/cristallina/analysis.py
+++ b/src/cristallina/analysis.py
@@ -326,6 +326,95 @@ def perform_image_roi_crop(
     return np.array(rois_within_batch)
 
 
+@memory.cache(ignore=["batch_size"])  # we ignore batch_size for caching purposes
+def calculate_image_histograms(
+    filesets,
+    channel="JF16T03V01",
+    alignment_channels=None,
+    batch_size=10,
+    roi: Optional[ROI] = None,
+    preview=False,
+    lower_cutoff_threshold=None,
+    bins=None,
+):
+    """
+    Calculates a histogram for a given region of interest (roi)
+    for an image channel from a fileset (e.g. "run0352/data/acq0001.*.h5" or step.fnames from a SFScanInfo object).
+
+    Allows alignment, i.e. reducing only to a common subset with other channels.
+
+    Calculations are performed in batches to reduce maximum memory requirements.
+
+    Preview only applies calculation to first batch and returns.
+
+    Returns:
+        (histogram, bins)
+    """
+
+    with SFDataFiles(*filesets) as data:
+        if alignment_channels is not None:
+            channels = [channel] + [ch for ch in alignment_channels]
+        else:
+            channels = [channel]
+
+        subset = data[channels]
+        subset.drop_missing()
+
+        Images = subset[channel]
+
+        # create empty array for stack sum with right shape
+        im = Images[0]
+        if roi is None:
+            im_ROI = im[:]
+        else:
+            im_ROI = im[roi.rows, roi.cols]
+
+        summed = np.zeros(im_ROI[0].shape)
+
+        if bins is None:
+            for image_slice in Images.in_batches(batch_size):
+                
+                index_slice, im = image_slice
+    
+                if roi is None:
+                    im_ROI = im
+                else:
+                    im_ROI = im[:, roi.rows, roi.cols]
+    
+                if lower_cutoff_threshold is not None:
+                    im_ROI = np.where(im_ROI < lower_cutoff_threshold, 0, im_ROI)
+                
+                bins = np.histogram_bin_edges(im.flatten(), bins='auto')
+                # only return first batch to calculate bins
+                break
+
+        all_hist = np.zeros(len(bins)-1)
+        for image_slice in Images.in_batches(batch_size):
+            
+            index_slice, im = image_slice
+
+            if roi is None:
+                im_ROI = im
+            else:
+                im_ROI = im[:, roi.rows, roi.cols]
+
+            if lower_cutoff_threshold is not None:
+                im_ROI = np.where(im_ROI < lower_cutoff_threshold, 0, im_ROI)
+            if bins is None:
+                bins = np.histogram_bin_edges(im.flatten(), bins='auto')
+            summed = summed + np.sum(im_ROI, axis=(0))
+
+            hist, _ = np.histogram(im.flatten(), bins=bins)
+            all_hist += hist
+            # only return first batch
+            if preview:
+                break
+    
+    return all_hist, bins
+
+
+
+
 def fit_2d_gaussian(image, roi: Optional[ROI] = None, plot=False):
     """
     2D Gaussian fit using LMFit for a given image and an optional region of interest.
@@ -585,4 +674,4 @@ def fit_1d_gaussian(x, y, use_offset=True, ax=None, print_results=False):
     if ax is not None:
         ax.plot(x, result.best_fit, label='fit')
 
-    return result
\ No newline at end of file
+    return result