diff --git a/shapeworks_cloud/core/migrations/0035_protect_from_cascading_deletion.py b/shapeworks_cloud/core/migrations/0035_protect_from_cascading_deletion.py
new file mode 100644
index 00000000..6f3da210
--- /dev/null
+++ b/shapeworks_cloud/core/migrations/0035_protect_from_cascading_deletion.py
@@ -0,0 +1,91 @@
+# Generated by Django 3.2.20 on 2023-09-10 04:58
+
+from django.conf import settings
+from django.db import migrations, models
+import django.db.models.deletion
+
+
+class Migration(migrations.Migration):
+    dependencies = [
+        migrations.swappable_dependency(settings.AUTH_USER_MODEL),
+        ('core', '0034_good_bad_particles'),
+    ]
+
+    operations = [
+        migrations.AlterField(
+            model_name='constraints',
+            name='optimized_particles',
+            field=models.ForeignKey(
+                null=True,
+                on_delete=django.db.models.deletion.SET_NULL,
+                related_name='constraints',
+                to='core.optimizedparticles',
+            ),
+        ),
+        migrations.AlterField(
+            model_name='dataset',
+            name='creator',
+            field=models.ForeignKey(
+                null=True, on_delete=django.db.models.deletion.SET_NULL, to=settings.AUTH_USER_MODEL
+            ),
+        ),
+        migrations.AlterField(
+            model_name='groomedmesh',
+            name='mesh',
+            field=models.ForeignKey(
+                null=True,
+                on_delete=django.db.models.deletion.SET_NULL,
+                related_name='groomed',
+                to='core.mesh',
+            ),
+        ),
+        migrations.AlterField(
+            model_name='groomedsegmentation',
+            name='segmentation',
+            field=models.ForeignKey(
+                null=True,
+                on_delete=django.db.models.deletion.SET_NULL,
+                related_name='groomed',
+                to='core.segmentation',
+            ),
+        ),
+        migrations.AlterField(
+            model_name='optimizedparticles',
+            name='groomed_mesh',
+            field=models.ForeignKey(
+                blank=True,
+                null=True,
+                on_delete=django.db.models.deletion.SET_NULL,
+                related_name='+',
+                to='core.groomedmesh',
+            ),
+        ),
+        migrations.AlterField(
+            model_name='optimizedparticles',
+            name='groomed_segmentation',
+            field=models.ForeignKey(
+                blank=True,
+                null=True,
+                on_delete=django.db.models.deletion.SET_NULL,
+                related_name='+',
+                to='core.groomedsegmentation',
+            ),
+        ),
+        migrations.AlterField(
+            model_name='project',
+            name='creator',
+            field=models.ForeignKey(
+                null=True, on_delete=django.db.models.deletion.SET_NULL, to=settings.AUTH_USER_MODEL
+            ),
+        ),
+        migrations.AlterField(
+            model_name='reconstructedsample',
+            name='particles',
+            field=models.ForeignKey(
+                null=True,
+                on_delete=django.db.models.deletion.SET_NULL,
+                related_name='reconstructed_samples',
+                to='core.optimizedparticles',
+            ),
+        ),
+    ]
diff --git a/shapeworks_cloud/core/migrations/0036_analysis_multi_domain.py b/shapeworks_cloud/core/migrations/0036_analysis_multi_domain.py
new file mode 100644
index 00000000..6c5f8f1a
--- /dev/null
+++ b/shapeworks_cloud/core/migrations/0036_analysis_multi_domain.py
@@ -0,0 +1,39 @@
+# Generated by Django 3.2.21 on 2023-09-25 20:42
+
+from django.db import migrations, models
+import s3_file_field.fields
+
+
+class Migration(migrations.Migration):
+    dependencies = [
+        ('core', '0035_protect_from_cascading_deletion'),
+    ]
+
+    operations = [
+        migrations.CreateModel(
+            name='CachedAnalysisMeanShape',
+            fields=[
+                (
+                    'id',
+                    models.AutoField(
+                        auto_created=True, primary_key=True, serialize=False, verbose_name='ID'
+                    ),
+                ),
+                ('file', s3_file_field.fields.S3FileField()),
+                ('particles', s3_file_field.fields.S3FileField(null=True)),
+            ],
+        ),
+        migrations.RemoveField(
+            model_name='cachedanalysis',
+            name='mean_particles',
+        ),
+        migrations.RemoveField(
+            model_name='cachedanalysis',
+            name='mean_shape',
+        ),
+        migrations.AddField(
+            model_name='cachedanalysis',
+            name='mean_shapes',
+            field=models.ManyToManyField(to='core.CachedAnalysisMeanShape'),
+        ),
+    ]
diff --git a/shapeworks_cloud/core/models.py b/shapeworks_cloud/core/models.py
index 6f638e3e..83210457 100644
--- a/shapeworks_cloud/core/models.py
+++ b/shapeworks_cloud/core/models.py
@@ -10,7 +10,7 @@
 class Dataset(TimeStampedModel, models.Model):
     name = models.CharField(max_length=255, unique=True)
     private = models.BooleanField(default=False)
-    creator = models.ForeignKey(User, on_delete=models.PROTECT, null=True)
+    creator = models.ForeignKey(User, on_delete=models.SET_NULL, null=True)
     thumbnail = S3FileField(null=True, blank=True)
     license = models.TextField()
     description = models.TextField()
@@ -95,9 +95,13 @@ class CachedAnalysisMode(models.Model):
     pca_values = models.ManyToManyField(CachedAnalysisModePCA)
 
 
+class CachedAnalysisMeanShape(models.Model):
+    file = S3FileField()
+    particles = S3FileField(null=True)
+
+
 class CachedAnalysis(TimeStampedModel, models.Model):
-    mean_shape = S3FileField()
-    mean_particles = S3FileField(null=True)
+    mean_shapes = models.ManyToManyField(CachedAnalysisMeanShape)
     modes = models.ManyToManyField(CachedAnalysisMode)
     charts = models.JSONField()
     groups = models.ManyToManyField(CachedAnalysisGroup, blank=True)
@@ -109,13 +113,17 @@ class Project(TimeStampedModel, models.Model):
     name = models.CharField(max_length=255)
     private = models.BooleanField(default=False)
     readonly = models.BooleanField(default=False)
-    creator = models.ForeignKey(User, on_delete=models.PROTECT, null=True)
+    creator = models.ForeignKey(User, on_delete=models.SET_NULL, null=True)
     thumbnail = S3FileField(null=True, blank=True)
     keywords = models.CharField(max_length=255, blank=True, default='')
     description = models.TextField(blank=True, default='')
     dataset = models.ForeignKey(Dataset, on_delete=models.CASCADE, related_name='projects')
-    last_cached_analysis = models.ForeignKey(CachedAnalysis, on_delete=models.SET_NULL, null=True)
     landmarks_info = models.JSONField(default=list, null=True)
+    last_cached_analysis = models.ForeignKey(
+        CachedAnalysis,
+        on_delete=models.SET_NULL,
+        null=True,
+    )
 
     def create_new_file(self):
         file_contents = {
@@ -151,15 +159,17 @@ def get_download_paths(self):
                     'groomed': [
                         (gm.mesh.anatomy_type, gm.file)
                         for gm in GroomedMesh.objects.filter(project=self, mesh__subject=subject)
+                        if gm.mesh
                     ]
                     + [
                         (gs.segmentation.anatomy_type, gs.file)
                         for gs in GroomedSegmentation.objects.filter(
                             project=self, segmentation__subject=subject
                         )
+                        if gs.segmentation
                     ],
-                    'local': [(p.anatomy_type, p.local) for p in particles],
-                    'world': [(p.anatomy_type, p.world) for p in particles],
+                    'local': [(p.anatomy_type, p.local) for p in particles if p.local],
+                    'world': [(p.anatomy_type, p.world) for p in particles if p.world],
                 }
                 related_files['shape'] = (
                     related_files['mesh']
@@ -176,7 +186,7 @@ def get_download_paths(self):
                         for related in related_files[prefix]:
                             if not target_file:
                                 # subject and anatomy type must match
-                                if suffix in related[0]:
+                                if suffix == related[0].replace('anatomy_', ''):
                                     target_file = related[1].url
                     if target_file:
                         value = value.replace('../', '')
@@ -194,8 +204,9 @@ class GroomedSegmentation(TimeStampedModel, models.Model):
 
     segmentation = models.ForeignKey(
         Segmentation,
-        on_delete=models.CASCADE,
+        on_delete=models.SET_NULL,
         related_name='groomed',
+        null=True,
     )
 
     project = models.ForeignKey(
@@ -213,8 +224,9 @@ class GroomedMesh(TimeStampedModel, models.Model):
 
     mesh = models.ForeignKey(
         Mesh,
-        on_delete=models.CASCADE,
+        on_delete=models.SET_NULL,
         related_name='groomed',
+        null=True,
     )
 
     project = models.ForeignKey(Project, on_delete=models.CASCADE, related_name='groomed_meshes')
@@ -232,14 +244,14 @@ class OptimizedParticles(TimeStampedModel, models.Model):
 
     groomed_segmentation = models.ForeignKey(
         GroomedSegmentation,
-        on_delete=models.CASCADE,
+        on_delete=models.SET_NULL,
         related_name='+',
         blank=True,
         null=True,
     )
     groomed_mesh = models.ForeignKey(
         GroomedMesh,
-        on_delete=models.CASCADE,
+        on_delete=models.SET_NULL,
         related_name='+',
         blank=True,
         null=True,
@@ -260,7 +272,10 @@ class Constraints(TimeStampedModel, models.Model):
     subject = models.ForeignKey(Subject, on_delete=models.CASCADE, related_name='constraints')
     anatomy_type = models.CharField(max_length=255)
     optimized_particles = models.ForeignKey(
-        OptimizedParticles, on_delete=models.CASCADE, related_name='constraints', null=True
+        OptimizedParticles,
+        on_delete=models.SET_NULL,
+        related_name='constraints',
+        null=True,
     )
 
 
@@ -270,7 +285,10 @@ class ReconstructedSample(TimeStampedModel, models.Model):
         Project, on_delete=models.CASCADE, related_name='reconstructed_samples'
     )
     particles = models.ForeignKey(
-        OptimizedParticles, on_delete=models.CASCADE, related_name='reconstructed_samples'
+        OptimizedParticles,
+        on_delete=models.SET_NULL,
+        related_name='reconstructed_samples',
+        null=True,
     )
 
 
diff --git a/shapeworks_cloud/core/rest.py b/shapeworks_cloud/core/rest.py
index 81dbf32d..7e6011fe 100644
--- a/shapeworks_cloud/core/rest.py
+++ b/shapeworks_cloud/core/rest.py
@@ -460,6 +460,11 @@ class CachedAnalysisGroupViewSet(BaseViewSet):
     serializer_class = serializers.CachedAnalysisGroupSerializer
 
 
+class CachedAnalysisMeanShapeViewSet(BaseViewSet):
+    queryset = models.CachedAnalysisMeanShape.objects.all()
+    serializer_class = serializers.CachedAnalysisMeanShapeSerializer
+
+
 class ReconstructedSampleViewSet(
     GenericViewSet,
     mixins.ListModelMixin,
diff --git a/shapeworks_cloud/core/serializers.py b/shapeworks_cloud/core/serializers.py
index cff33cd1..267ee33c 100644
--- a/shapeworks_cloud/core/serializers.py
+++ b/shapeworks_cloud/core/serializers.py
@@ -18,6 +18,12 @@ class Meta:
         fields = '__all__'
 
 
+class CachedAnalysisMeanShapeSerializer(serializers.ModelSerializer):
+    class Meta:
+        model = models.CachedAnalysisMeanShape
+        fields = '__all__'
+
+
 class CachedAnalysisModePCASerializer(serializers.ModelSerializer):
     class Meta:
         model = models.CachedAnalysisModePCA
@@ -55,6 +61,7 @@ class Meta:
 class CachedAnalysisReadSerializer(serializers.ModelSerializer):
     modes = CachedAnalysisModeReadSerializer(many=True)
     groups = CachedAnalysisGroupSerializer(many=True)
+    mean_shapes = CachedAnalysisMeanShapeSerializer(many=True)
 
     class Meta:
         model = models.CachedAnalysis
@@ -112,6 +119,18 @@ class Meta:
 
 
 class SubjectSerializer(serializers.ModelSerializer):
+    num_domains = serializers.SerializerMethodField('get_num_domains')
+
+    def get_num_domains(self, obj):
+        shapes = list(obj.segmentations.all()) + list(obj.meshes.all()) + list(obj.contours.all())
+        domains = []
+        for shape in shapes:
+            # get unique values for anatomy_type
+            domain = shape.anatomy_type
+            if domain not in domains:
+                domains.append(domain)
+        return len(domains)
+
     class Meta:
         model = models.Subject
         fields = '__all__'
diff --git a/shapeworks_cloud/core/signals.py b/shapeworks_cloud/core/signals.py
index 27e6259c..72e43357 100644
--- a/shapeworks_cloud/core/signals.py
+++ b/shapeworks_cloud/core/signals.py
@@ -1,7 +1,14 @@
 from django.db.models.signals import pre_delete
 from django.dispatch import receiver
 
-from .models import CachedAnalysis, CachedAnalysisMode, CachedAnalysisModePCA, Project
+from .models import (
+    CachedAnalysis,
+    CachedAnalysisGroup,
+    CachedAnalysisMeanShape,
+    CachedAnalysisMode,
+    CachedAnalysisModePCA,
+    Project,
+)
 
 
 @receiver(pre_delete, sender=Project)
@@ -11,3 +18,5 @@ def delete_cached_analysis(sender, instance, using, **kwargs):
     ).delete()
     CachedAnalysisMode.objects.filter(cachedanalysis__project=instance).delete()
     CachedAnalysis.objects.filter(project=instance).delete()
+    CachedAnalysisGroup.objects.filter(cachedanalysis__project=instance).delete()
+    CachedAnalysisMeanShape.objects.filter(cachedanalysis__project=instance).delete()
diff --git a/shapeworks_cloud/core/tasks.py b/shapeworks_cloud/core/tasks.py
index e800290a..bc7ab183 100644
--- a/shapeworks_cloud/core/tasks.py
+++ b/shapeworks_cloud/core/tasks.py
@@ -10,7 +10,6 @@
 from django.contrib.auth.models import User
 from django.core.files.base import ContentFile
 from django.db.models import Q
-import pandas
 from rest_framework.authtoken.models import Token
 
 from shapeworks_cloud.core import models
@@ -28,10 +27,9 @@ def parse_progress(xml_string):
     return 0
 
 
-def edit_swproj_section(filename, section_name, new_df):
+def edit_swproj_section(filename, section_name, new_contents):
     with open(filename, 'r') as f:
         data = json.load(f)
-    new_contents = {item['key']: item['value'] for item in new_df.to_dict(orient='records')}
     if section_name == 'groom':
         data[section_name] = {}
         data[section_name]['shape'] = new_contents
@@ -48,28 +46,28 @@ def edit_swproj_section(filename, section_name, new_df):
         json.dump(data, f)
 
 
-def interpret_form_df(df, command):
-    if command == 'groom' and df['key'].str.contains('anisotropic_').any():
-        # consolidate anisotropic values to one row
-        anisotropic_values = {
-            axis: str(df.loc[df['key'] == 'anisotropic_' + axis].iloc[0]['value'])
-            for axis in ['x', 'y', 'z']
-        }
-        df_filter = df['key'].map(lambda key: 'anisotropic_' not in key)
-        df = df[df_filter]
-        return pandas.concat(
-            [
-                df,
-                pandas.DataFrame.from_dict(
-                    {
-                        'key': ['spacing'],
-                        'value': [' '.join(anisotropic_values.values())],
-                    }
-                ),
-            ]
-        )
-    else:
-        return df
+def interpret_form_data(data, command, swcc_project):
+    anisotropic_values = []
+    del_keys = []
+    for key, value in data.items():
+        if 'anisotropic' in key:
+            anisotropic_values.append(value)
+            del_keys.append(key)
+
+    for del_key in del_keys:
+        del data[del_key]
+
+    if command == 'groom' and len(anisotropic_values) > 0:
+        data['spacing'] = ' '.join(anisotropic_values)
+    elif command == 'optimize':
+        num_particles = data.get('number_of_particles')
+        if num_particles:
+            max_num_domains = max(s.num_domains for s in swcc_project.subjects)
+            data['number_of_particles'] = ' '.join(
+                str(num_particles) for i in range(max_num_domains)
+            )
+
+    return data
 
 
 def run_shapeworks_command(
@@ -93,22 +91,19 @@ def run_shapeworks_command(
             session.set_token(token.key)
             project = models.Project.objects.get(id=project_id)
             project_filename = project.file.name.split('/')[-1]
-            SWCCProject.from_id(project.id).download(download_dir)
+            swcc_project = SWCCProject.from_id(project.id)
+            swcc_project.download(download_dir)
 
             pre_command_function()
             progress.update_percentage(10)
 
             if form_data:
                 # write the form data to the project file
-                form_df = pandas.DataFrame(
-                    list(form_data.items()),
-                    columns=['key', 'value'],
-                )
-                form_df = interpret_form_df(form_df, command)
+                form_data = interpret_form_data(form_data, command, swcc_project)
                 edit_swproj_section(
                     Path(download_dir, project_filename),
                     command,
-                    form_df,
+                    form_data,
                 )
 
             # perform command
@@ -177,12 +172,9 @@ def post_command_function(project, download_dir, result_data, project_filename):
                 if len(prefixes) > 0:
                     prefix = prefixes[0]
                     anatomy_id = 'anatomy' + key.replace(prefix, '')
-                    if prefix in ['mesh', 'segmentation', 'image', 'contour']:
-                        prefix = 'shape'
-                    if prefix in ['shape', 'groomed']:
-                        if anatomy_id not in row:
-                            row[anatomy_id] = {}
-                        row[anatomy_id][prefix] = entry[key].replace('../', '').replace('./', '')
+                    if anatomy_id not in row:
+                        row[anatomy_id] = {}
+                    row[anatomy_id][prefix] = entry[key].replace('../', '').replace('./', '')
 
             for anatomy_data in row.values():
                 if 'groomed' not in anatomy_data:
@@ -262,9 +254,9 @@ def post_command_function(project, download_dir, result_data, project_filename):
                 target_mesh = project_groomed_meshes.filter(
                     file__endswith=groomed_filename,
                 ).first()
-                if target_mesh:
+                if target_mesh and target_mesh.mesh:
                     subject = target_mesh.mesh.subject
-                elif target_segmentation:
+                elif target_segmentation and target_segmentation.segmentation:
                     subject = target_segmentation.segmentation.subject
                 result_particles_object = models.OptimizedParticles.objects.create(
                     groomed_segmentation=target_segmentation,
@@ -314,6 +306,8 @@ def pre_command_function():
             cachedanalysismode__cachedanalysis__project=project
         ).delete()
         models.CachedAnalysisMode.objects.filter(cachedanalysis__project=project).delete()
+        models.CachedAnalysisGroup.objects.filter(cachedanalysis__project=project).delete()
+        models.CachedAnalysisMeanShape.objects.filter(cachedanalysis__project=project).delete()
         models.CachedAnalysis.objects.filter(project=project).delete()
 
     def post_command_function(project, download_dir, result_data, project_filename):
@@ -328,16 +322,16 @@ def post_command_function(project, download_dir, result_data, project_filename):
             project_data = json.load(pf)['data']
             for i, sample in enumerate(project_data):
                 reconstructed_filenames = result_data['reconstructed_samples'][i]
-                particles = (
-                    models.OptimizedParticles.objects.filter(project=project)
-                    .filter(
+                subject_particles = list(
+                    models.OptimizedParticles.objects.filter(project=project).filter(
                         Q(groomed_mesh__mesh__subject__name=sample['name'])
                         | Q(groomed_segmentation__segmentation__subject__name=sample['name'])
                     )
-                    .first()
                 )
-                for reconstructed_filename in reconstructed_filenames:
-                    reconstructed = models.ReconstructedSample(project=project, particles=particles)
+                for j, reconstructed_filename in enumerate(reconstructed_filenames):
+                    reconstructed = models.ReconstructedSample(
+                        project=project, particles=subject_particles[j]
+                    )
                     reconstructed.file.save(
                         reconstructed_filename,
                         open(Path(download_dir, reconstructed_filename), 'rb'),
diff --git a/shapeworks_cloud/urls.py b/shapeworks_cloud/urls.py
index 51344c9a..23daeb07 100644
--- a/shapeworks_cloud/urls.py
+++ b/shapeworks_cloud/urls.py
@@ -41,6 +41,11 @@
     rest.CachedAnalysisGroupViewSet,
     basename='cached_analysis_group',
 )
+router.register(
+    'cached-analysis-mean-shape',
+    rest.CachedAnalysisMeanShapeViewSet,
+    basename='cached_analysis_mean_shape',
+)
 router.register(
     'reconstructed-samples', rest.ReconstructedSampleViewSet, basename='reconstructed_sample'
 )
diff --git a/swcc/swcc/models/constants.py b/swcc/swcc/models/constants.py
index d4b739be..5a98f1e5 100644
--- a/swcc/swcc/models/constants.py
+++ b/swcc/swcc/models/constants.py
@@ -1,5 +1,12 @@
+required_key_prefixes = [
+    'shape',
+    'mesh',
+    'segmentation',
+    'contour',
+    'image',
+]
+
 expected_key_prefixes = [
-    'name',
     'shape',
     'mesh',
     'segmentation',
@@ -8,8 +15,5 @@
     'groomed',
     'local',
     'world',
-    'alignment',
-    'procrustes',
-    'landmarks',
     'constraints',
 ]
diff --git a/swcc/swcc/models/other_models.py b/swcc/swcc/models/other_models.py
index b6b3ab95..f35d9f50 100644
--- a/swcc/swcc/models/other_models.py
+++ b/swcc/swcc/models/other_models.py
@@ -131,11 +131,17 @@ class CachedAnalysisMode(ApiModel):
     pca_values: List[CachedAnalysisModePCA]
 
 
+class CachedAnalysisMeanShape(ApiModel):
+    _endpoint = 'cached-analysis-mean-shape'
+
+    file: FileType[Literal['core.CachedAnalysisMeanShape.file']]
+    particles: FileType[Literal['core.CachedAnalysisMeanShape.particles']]
+
+
 class CachedAnalysis(ApiModel):
     _endpoint = 'cached-analysis'
 
-    mean_shape: FileType[Literal['core.CachedAnalysis.mean_shape']]
-    mean_particles: FileType[Literal['core.CachedAnalysis.mean_particles']]
+    mean_shapes: List[CachedAnalysisMeanShape]
     modes: List[CachedAnalysisMode]
     charts: List[dict]
     groups: Optional[List[CachedAnalysisGroup]]
diff --git a/swcc/swcc/models/project.py b/swcc/swcc/models/project.py
index b8fa7efa..2ae204a8 100644
--- a/swcc/swcc/models/project.py
+++ b/swcc/swcc/models/project.py
@@ -3,6 +3,7 @@
 import json
 from pathlib import Path
 from tempfile import TemporaryDirectory
+import warnings
 
 import requests
 
@@ -24,12 +25,13 @@
 
 from ..api import current_session
 from .api_model import ApiModel
-from .constants import expected_key_prefixes
+from .constants import expected_key_prefixes, required_key_prefixes
 from .dataset import Dataset
 from .file_type import FileType
 from .other_models import (
     CachedAnalysis,
     CachedAnalysisGroup,
+    CachedAnalysisMeanShape,
     CachedAnalysisMode,
     CachedAnalysisModePCA,
     Constraints,
@@ -97,25 +99,26 @@ def interpret_data(self, input_data):
                     name=entry.get('name'), groups=groups_dict, dataset=self.project.dataset
                 ).create()
 
-            entry_values: Dict = {p: [] for p in expected_key_prefixes}
-            entry_values['anatomy_ids'] = []
+            objects_by_domain: Dict[str, Dict] = {}
             for key in entry.keys():
-                if key != 'name':
-                    prefixes = [p for p in expected_key_prefixes if key.startswith(p)]
-                    if len(prefixes) > 0:
-                        entry_values[prefixes[0]].append(entry[key])
-                        anatomy_id = 'anatomy' + key.replace(prefixes[0], '').replace(
-                            '_particles', ''
-                        ).replace('_file', '')
-                        if anatomy_id not in entry_values['anatomy_ids']:
-                            entry_values['anatomy_ids'].append(anatomy_id)
-            objects_by_domain = {}
-            for index, anatomy_id in enumerate(entry_values['anatomy_ids']):
-                objects_by_domain[anatomy_id] = {
-                    k: v[index] if len(v) > index else v[0]
-                    for k, v in entry_values.items()
-                    if len(v) > 0
-                }
+                prefixes = [p for p in expected_key_prefixes if key.startswith(p)]
+                if len(prefixes) > 0:
+                    prefix = prefixes[0]
+                    anatomy_id = 'anatomy' + key
+                    anatomy_id = anatomy_id.replace(prefix, '').replace('_particles', '')
+                    # Only create a new domain object if a shape exists for that suffix
+                    if anatomy_id not in objects_by_domain:
+                        if prefix in required_key_prefixes:
+                            objects_by_domain[anatomy_id] = {}
+                        else:
+                            warnings.warn(
+                                f'No shape exists for {anatomy_id}. Cannot create {key}.',
+                                stacklevel=2,
+                            )
+                            continue
+                    objects_by_domain[anatomy_id][prefix] = (
+                        entry[key].replace('../', '').replace('./', '')
+                    )
             output_data.append(
                 [
                     subject,
@@ -235,12 +238,23 @@ def load_analysis_from_json(self, file_path):
         analysis_file_location = project_root / Path(file_path)
         contents = json.load(open(analysis_file_location))
         if contents['mean'] and contents['mean']['meshes']:
-            mean_shape_path = contents['mean']['meshes'][0]
-            mean_particles_path = None
+            mean_shapes_cache = []
+            mean_shapes = []
+            for mean_shape in contents['mean']['meshes']:
+                mean_shapes.append(analysis_file_location.parent / Path(mean_shape))
+
             if 'particle_files' in contents['mean']:
-                mean_particles_path = contents['mean']['particle_files'][0]
-            if 'particles' in contents['mean']:
-                mean_particles_path = contents['mean']['particles'][0]
+                mean_particles = []
+                for mean_particle_path in contents['mean']['particle_files']:
+                    mean_particles.append(analysis_file_location.parent / Path(mean_particle_path))
+
+            for i in range(len(mean_shapes)):
+                cams = CachedAnalysisMeanShape(
+                    file=mean_shapes[i],
+                    particles=mean_particles[i] if mean_particles else None,
+                ).create()
+                mean_shapes_cache.append(cams)
+
             modes = []
             for mode in contents['modes']:
                 pca_values = []
@@ -271,10 +285,6 @@ def load_analysis_from_json(self, file_path):
                     modes.append(cam)
 
             if len(modes) > 0:
-                mean_particles = None
-                if mean_particles_path:
-                    mean_particles = analysis_file_location.parent / Path(mean_particles_path)
-
                 groups_cache = []
                 if contents['groups']:
                     for group in contents['groups']:
@@ -295,8 +305,7 @@ def load_analysis_from_json(self, file_path):
 
                                 groups_cache.append(cag)
                 return CachedAnalysis(
-                    mean_shape=analysis_file_location.parent / Path(mean_shape_path),
-                    mean_particles=mean_particles,
+                    mean_shapes=mean_shapes_cache,
                     modes=modes,
                     charts=contents['charts'],
                     groups=groups_cache,
@@ -322,6 +331,10 @@ class Project(ApiModel):
     def get_file_io(self):
         return ProjectFileIO(project=self)
 
+    @property
+    def subjects(self) -> Iterator[Subject]:
+        return Subject.list(project=self)
+
     @property
     def groomed_segmentations(self) -> Iterator[GroomedSegmentation]:
         self.assert_remote()
diff --git a/swcc/swcc/models/subject.py b/swcc/swcc/models/subject.py
index 1fc3a18f..334ef3dd 100644
--- a/swcc/swcc/models/subject.py
+++ b/swcc/swcc/models/subject.py
@@ -13,6 +13,7 @@ class Subject(ApiModel):
     name: NonEmptyString
     dataset: Dataset
     groups: Optional[Dict[str, str]]
+    num_domains: Optional[int]
 
     @property
     def segmentations(self) -> Iterator[Segmentation]:
diff --git a/swcc/tests/test_download_upload.py b/swcc/tests/test_download_upload.py
index aeb9d8b3..faff1b43 100644
--- a/swcc/tests/test_download_upload.py
+++ b/swcc/tests/test_download_upload.py
@@ -49,12 +49,12 @@ def is_same(dir1, dir2):
 def public_server_download(download_dir):
     with swcc_session() as public_server_session:
         public_server_session.login('testuser@noemail.nil', 'cicdtest')
-        all_projects = list(models.Project.list())
-        project_subset = (
-            random.sample(all_projects, SAMPLE_SIZE)
-            if len(all_projects) >= SAMPLE_SIZE
-            else all_projects
+        all_datasets = list(models.Dataset.list())
+        tiny_tests = [d for d in all_datasets if 'tiny_test' in d.name]
+        dataset_subset = (
+            random.sample(tiny_tests, SAMPLE_SIZE) if len(tiny_tests) >= SAMPLE_SIZE else tiny_tests
         )
+        project_subset = [next(d.projects) for d in dataset_subset]
         for project in project_subset:
             project.download(download_dir)
         return project_subset
diff --git a/web/shapeworks/src/components/Analysis/AnalysisTab.vue b/web/shapeworks/src/components/Analysis/AnalysisTab.vue
index c1cc0913..e1d764bd 100644
--- a/web/shapeworks/src/components/Analysis/AnalysisTab.vue
+++ b/web/shapeworks/src/components/Analysis/AnalysisTab.vue
@@ -1,7 +1,6 @@
 <script lang="ts">
 import {
     analysis,
-    analysisFileShown,
     currentTasks,
     selectedProject,
     analysisExpandedTab,
@@ -79,7 +78,6 @@ export default {
 
         return {
             analysis,
-            analysisFileShown,
             message,
             taskData,
             openTab,
diff --git a/web/shapeworks/src/components/Analysis/Groups.vue b/web/shapeworks/src/components/Analysis/Groups.vue
index 428d7d57..90b1e9f4 100644
--- a/web/shapeworks/src/components/Analysis/Groups.vue
+++ b/web/shapeworks/src/components/Analysis/Groups.vue
@@ -1,7 +1,13 @@
 <script lang="ts">
-  import { analysis, analysisFileShown, cacheAllComparisons, currentAnalysisFileParticles, meanAnalysisFileParticles } from '@/store';
+import {
+    analysis,
+    cacheAllComparisons,
+    analysisFilesShown,
+    currentAnalysisParticlesFiles,
+    meanAnalysisParticlesFiles
+} from '@/store';
 import { AnalysisGroup } from '@/types';
-import { Ref, computed, inject, ref, watch } from 'vue';
+import { Ref, computed, inject, ref, watch, onMounted } from 'vue';
 
 import { AnalysisTabs } from './util';
 
@@ -11,7 +17,7 @@ import { AnalysisTabs } from './util';
       openTab: Number,
     },
     setup(props) {
-      const currGroup = ref<AnalysisGroup>();
+      const currGroup = ref<AnalysisGroup[]>();
       const groupRatio = ref(0.5);
       const groupDiff = ref(false);
       const groupSet = ref<string>();
@@ -31,7 +37,7 @@ import { AnalysisTabs } from './util';
             return (g.group1 === currPairing.value.right && g.group2 === currPairing.value.left)
         },
         updateCurrGroup() {
-            currGroup.value = analysis.value?.groups.find((g) => {
+            currGroup.value = analysis.value?.groups.filter((g) => {
                 if (g.name === groupSet.value) {
                     if ((g.group1 === currPairing.value.left && g.group2 === currPairing.value.right)
                         || (g.group1 === currPairing.value.right && g.group2 === currPairing.value.left)) { // if the groups pairings match the selected pairings
@@ -49,6 +55,9 @@ import { AnalysisTabs } from './util';
             })
         },
         setDefaultPairing() {
+            if (groupSet.value === undefined && allGroupSets.value?.length) {
+                groupSet.value = allGroupSets.value[0];
+            }
             // default left and right group selections: first and second item in groupSet pairings list
             if (groupSet.value !== undefined) {
                 currPairing.value.left = groupPairings.value[groupSet.value][0];
@@ -75,17 +84,17 @@ import { AnalysisTabs } from './util';
             methods.updateGroupFileShown();
         },
         updateGroupFileShown() {
-            let fileShown: string | undefined = undefined;
-            let particles: string | undefined = undefined;
+            let filesShown: string[] | undefined = undefined;
+            let particles: string[] | undefined = undefined;
             methods.updateCurrGroup();
 
             if (props.currentTab === 'analyze' && analysis.value && currGroup.value){
-                fileShown = currGroup.value.file;
-                particles = currGroup.value.particles;
+                filesShown = currGroup.value.map((g) => g.file);
+                particles = currGroup.value.map((g) => g.particles);
             }
-            currentAnalysisFileParticles.value = particles;
+            currentAnalysisParticlesFiles.value = particles;
             if (groupDiff.value) {
-                meanAnalysisFileParticles.value = analysis.value?.groups.find((g) => {
+                const meanParticles = analysis.value?.groups.filter((g) => {
                     if (g.name === groupSet.value) {
                       if (!methods.pairingIsInverted(g)) {
                         if (g.ratio === 1.0) return true;
@@ -93,12 +102,15 @@ import { AnalysisTabs } from './util';
                         if (g.ratio === 0.0) return true;
                       }
                     }
-                })?.particles // account for inverted-pairings
+                }) // account for inverted-pairings
+                if (meanParticles) {
+                    meanAnalysisParticlesFiles.value = meanParticles.map((p) => p.particles)
+                }
             } else {
-                meanAnalysisFileParticles.value = analysis.value?.mean_particles;
+                meanAnalysisParticlesFiles.value = analysis.value?.mean_shapes.map((m) => m.particles)
             }
 
-            analysisFileShown.value = fileShown;
+            analysisFilesShown.value = filesShown;
         },
         async cacheAllGroupComparisons() {
             const allInPairing = analysis.value?.groups.filter((g) => {
@@ -109,7 +121,7 @@ import { AnalysisTabs } from './util';
             })
 
             if (allInPairing !== undefined) {
-                await cacheAllComparisons(allInPairing);
+                await cacheAllComparisons([allInPairing]);
             }
         },
         animateSlider() {
@@ -152,14 +164,15 @@ import { AnalysisTabs } from './util';
             return g;
         })
 
+        onMounted(methods.setDefaultPairing);
         watch(groupSet, () => {
             methods.setDefaultPairing();
             methods.updateGroupFileShown();
-            animate.value = false;
+            methods.stopAnimating();
         })
         watch(currPairing.value, () => {
             methods.updateGroupSelections();
-            animate.value = false;
+            methods.stopAnimating();
         })
         watch(groupRatio, methods.updateGroupFileShown)
         watch(groupDiff, methods.updateGroupFileShown)
@@ -178,7 +191,7 @@ import { AnalysisTabs } from './util';
         currGroup,
         currPairing,
         animate,
-        currentlyCaching
+        currentlyCaching,
       };
     },
   };
@@ -237,9 +250,9 @@ import { AnalysisTabs } from './util';
           ></v-checkbox>
       </v-row>
       <v-card align="center" justify="center" class="ma-auto mb-3" :disabled="currGroup === undefined">
-          <v-btn class="ms-4" color="grey darken-3" @click="() => { groupRatio = 0.0;}">Left Mean</v-btn>
+          <v-btn class="ms-4" color="grey darken-3" @click="() => { groupRatio = 0.0; methods.stopAnimating()}">Left Mean</v-btn>
           <v-btn-toggle class="ms-4" color="white"><v-btn color="grey darken-4" :disabled="animate || currentlyCaching" @click="() => groupDiff = !groupDiff">Diff --></v-btn></v-btn-toggle>
-          <v-btn class="ms-4" color="grey darken-3" @click="() => { groupRatio = 1.0;}">Right Mean</v-btn>
+          <v-btn class="ms-4" color="grey darken-3" @click="() => { groupRatio = 1.0; methods.stopAnimating()}">Right Mean</v-btn>
       </v-card>
     </div>
 </template>
diff --git a/web/shapeworks/src/components/Analysis/PCA.vue b/web/shapeworks/src/components/Analysis/PCA.vue
index a521c760..d0fa3774 100644
--- a/web/shapeworks/src/components/Analysis/PCA.vue
+++ b/web/shapeworks/src/components/Analysis/PCA.vue
@@ -1,14 +1,15 @@
 <script lang="ts">
   import { computed, ref, Ref, watch, inject } from 'vue';
   import { AnalysisTabs } from './util';
-  import { AnalysisParams } from '@/types';
+  import { AnalysisParams, CacheComparison } from '@/types';
+import { groupBy } from '../../helper';
   import {
     analysis,
-    analysisFileShown,
+    analysisFilesShown,
     cacheAllComparisons,
-    currentAnalysisFileParticles,
+    currentAnalysisParticlesFiles,
     currentTasks,
-    meanAnalysisFileParticles,
+    meanAnalysisParticlesFiles,
     selectedProject,
     spawnJob,
     spawnJobProgressPoll
@@ -42,8 +43,8 @@
         return analysis.value?.modes.find((m) => m.mode == mode.value)
       })
 
-      const currPCA = computed(() => {
-        return currMode.value?.pca_values.find(
+      const currPCAs = computed(() => {
+        return currMode.value?.pca_values.filter(
             p => p.pca_value.toPrecision(2) === stdDev.value.toPrecision(2)
         )
       })
@@ -51,10 +52,11 @@
       const stdDevRange = computed(() => {
         if (!analysis.value || !currMode.value) return {range: 0.0, step: 0.0, numSteps: 0}
         const pcaValues = currMode.value.pca_values.map(p => Math.round(p.pca_value * 100) / 100)
+        const numDomains = analysis.value.mean_shapes.length || 1
         const min =  Math.min(...pcaValues);
         const max =  Math.max(...pcaValues);
         const range = max;
-        const step = (max - min) / pcaValues.length;
+        const step = (max - min) / (pcaValues.length / numDomains);
         const numSteps = Math.round((max - min) / step) + 1;
 
         return {
@@ -71,7 +73,7 @@
                 {text: 'value', value: 'value', align: 'end'}
             ],
             items: [
-                {key: 'Lambda', value:currPCA.value ? Math.round(currPCA.value.lambda_value * 100) / 100 : 0},
+                {key: 'Lambda', value: currPCAs.value?.length ? Math.round(currPCAs.value[0].lambda_value * 100) / 100 : 0},
                 {key: 'Eigenvalue', value: currMode.value ? Math.round(currMode.value.eigen_value * 100) / 100 : 0},
                 {key: 'Explained Variance', value: currMode.value ? Math.round(currMode.value.explained_variance * 100) / 100 : 0, class: 'percentage'},
                 {key: 'Cumulative Explained Variance', value: currMode.value ? Math.round(currMode.value.cumulative_explained_variance * 100) / 100 : 0, class:'percentage'},
@@ -81,25 +83,24 @@
 
       const methods = {
         updateFileShown() {
-            let fileShown: string | undefined = undefined
-            let particles: string | undefined = undefined
+            let filesShown: string[] = []
+            let particles: string[] = []
             if (props.currentTab === 'analyze' && analysis.value){
-                if (stdDev.value === 0) {
-                    fileShown = analysis.value.mean_shape
-                    particles = analysis.value.mean_particles;
-                } else {
-                    fileShown = currPCA.value?.file
-                    particles = currPCA.value?.particles
-                }
+              if (stdDev.value === 0) {
+                filesShown = analysis.value.mean_shapes.map((m) => m.file)
+                particles = analysis.value.mean_shapes.map((m) => m.particles);
+              } else if (currPCAs.value) {
+                filesShown = currPCAs.value.map(p => p.file)
+                particles = currPCAs.value.map(p => p.particles)
+              }
             }
-            currentAnalysisFileParticles.value = particles;
-            meanAnalysisFileParticles.value =  analysis.value?.mean_particles;
-            analysisFileShown.value = fileShown;
+            currentAnalysisParticlesFiles.value = particles;
+            meanAnalysisParticlesFiles.value =  analysis.value?.mean_shapes.map((m) => m.particles);
+            analysisFilesShown.value = filesShown;
         },
         async cacheAllPCAComparisons() {
           if (currMode.value !== undefined) {
-            const allInMode = currMode.value.pca_values;
-
+            const allInMode: CacheComparison[][] = Object.values(groupBy(currMode.value.pca_values, 'pca_value'));
             await cacheAllComparisons(allInMode);
           }
         },
diff --git a/web/shapeworks/src/components/RenderControls.vue b/web/shapeworks/src/components/RenderControls.vue
index c8567f7f..5de56a9a 100644
--- a/web/shapeworks/src/components/RenderControls.vue
+++ b/web/shapeworks/src/components/RenderControls.vue
@@ -13,7 +13,7 @@ import {
     allDatasets,
     allProjectsForDataset,
     showDifferenceFromMeanMode,
-    analysisFileShown,
+    analysisFilesShown,
     showGoodBadParticlesMode,
 goodBadMaxAngle
 } from '@/store';
@@ -158,7 +158,7 @@ export default {
         }
 
         const showAnalysisOptions = computed(() => {
-            return props.currentTab === 'analyze' && analysisFileShown.value;
+            return props.currentTab === 'analyze' && analysisFilesShown.value?.length;
         })
 
         return {
diff --git a/web/shapeworks/src/components/ShapeViewer/methods.js b/web/shapeworks/src/components/ShapeViewer/methods.js
index 76ad5964..9f35d1dd 100644
--- a/web/shapeworks/src/components/ShapeViewer/methods.js
+++ b/web/shapeworks/src/components/ShapeViewer/methods.js
@@ -18,7 +18,7 @@ import vtkPolyData from 'vtk.js/Sources/Common/DataModel/PolyData';
 import {
     layers, layersShown, orientationIndicator,
     cachedMarchingCubes, cachedParticleComparisonColors, vtkShapesByType,
-    analysisFileShown, currentAnalysisFileParticles, meanAnalysisFileParticles,
+    analysisFilesShown, currentAnalysisParticlesFiles, meanAnalysisParticlesFiles,
     showDifferenceFromMeanMode, cachedParticleComparisonVectors, landmarkColorList, cacheComparison, calculateComparisons,
     showGoodBadParticlesMode, goodBadMaxAngle, goodBadAngles,
 } from '@/store';
@@ -117,7 +117,7 @@ export default {
             this.applyCameraDelta(renderer, positionDelta, viewUpDelta)
         })
     },
-    createColorFilter(landmarks = false, goodBad = false) {
+    createColorFilter(domainIndex = 0, landmarks = false, goodBad = false) {
         const filter = vtkCalculator.newInstance()
         filter.setFormula({
             getArrays() {
@@ -147,10 +147,8 @@ export default {
                         }
                     }
 
-                    if (goodBad && analysisFileShown.value) {
-                        // TODO: determine domain index for each particle array
-                        // MULTIDOMAIN FIX
-                        if (goodBadAngles.value[0][i] < goodBadMaxAngle.value) {
+                    if (goodBad && analysisFilesShown.value?.length) {
+                        if (goodBadAngles.value[domainIndex][i] < goodBadMaxAngle.value) {
                             c = GOOD_BAD_COLORS[0] // green
                         } else {
                             c = GOOD_BAD_COLORS[1] // red
@@ -165,7 +163,7 @@ export default {
         });
         return filter;
     },
-    addPoints(renderer, points, landmarks = false) {
+    addPoints(renderer, points, i, landmarks = false) {
         let size = this.glyphSize
         let source = vtkSphereSource.newInstance({
             thetaResolution: SPHERE_RESOLUTION,
@@ -182,7 +180,7 @@ export default {
             scaleFactor: size,
         });
         const actor = vtkActor.newInstance();
-        const filter = this.createColorFilter(landmarks, showGoodBadParticlesMode.value);
+        const filter = this.createColorFilter(i, landmarks, showGoodBadParticlesMode.value);
 
         filter.setInputData(points, 0);
         mapper.setInputConnection(filter.getOutputPort(), 0);
@@ -193,65 +191,71 @@ export default {
         this.vtk.pointMappers.push(mapper);
     },
     addShapes(renderer, label, shapes) {
-        shapes.flat().forEach(
-            (shape, index) => {
-                let layerName = Object.entries(vtkShapesByType.value).filter(
-                    ([, shapes]) => shapes.includes(shape)
-                ).map(
-                    ([layerName,]) => layerName
-                )
-                layerName = layerName.length ? layerName[0] : "Original"
-                const type = layers.value.find((layer) => layer.name === layerName)
-                let opacity = 1;
-                if (!analysisFileShown.value) {
-                    const numLayers = layersShown.value.filter(
-                        (layerName) => layers.value.find((layer) => layer.name == layerName).rgb
-                    ).length
-                    if (numLayers > 0) opacity /= numLayers
-                }
-                const cacheLabel = `${label}_${layerName}_${index}`
+        shapes.forEach(
+            (shapeDatas, domainIndex) => {
+                shapeDatas.forEach(
+                    (shapeData) => {
+                        let layerName = Object.entries(vtkShapesByType.value).filter(
+                            ([, shapes]) => shapes.includes(shapeData)
+                        ).map(
+                            ([layerName,]) => layerName
+                        )
+                        layerName = layerName.length ? layerName[0] : "Original"
+                        const type = layers.value.find((layer) => layer.name === layerName)
+                        let opacity = 1;
+                        if (!analysisFilesShown.value?.length) {
+                            const numLayers = layersShown.value.filter(
+                                (layerName) => layers.value.find((layer) => layer.name == layerName).rgb
+                            ).length
+                            if (numLayers > 0) opacity /= numLayers
+                        }
+                        const cacheLabel = `${label}_${layerName}_${domainIndex}`
 
-                const mapper = vtkMapper.newInstance({
-                    colorMode: ColorMode.MAP_SCALARS,
-                    scalarMode: ScalarMode.USE_POINT_FIELD_DATA,
-                });
-                const actor = vtkActor.newInstance();
-                actor.getProperty().setColor(...type.rgb);
-                actor.getProperty().setOpacity(opacity);
-                actor.setMapper(mapper);
-                if (shape.getClassName() == 'vtkPolyData') {
-                    mapper.setInputData(shape);
-                } else if (cachedMarchingCubes.value[cacheLabel]) {
-                    mapper.setInputData(cachedMarchingCubes.value[cacheLabel])
-                } else {
-                    const marchingCube = vtkImageMarchingCubes.newInstance({
-                        contourValue: 0.001,
-                        computeNormals: true,
-                        mergePoints: true,
-                    });
-                    marchingCube.setInputData(shape)
-                    mapper.setInputConnection(marchingCube.getOutputPort());
-                    cachedMarchingCubes.value[cacheLabel] = marchingCube.getOutputData()
-                }
-                if (showDifferenceFromMeanMode.value) {
-                    this.showDifferenceFromMean(mapper, renderer)
-                }
-                renderer.addActor(actor);
+                        const mapper = vtkMapper.newInstance({
+                            colorMode: ColorMode.MAP_SCALARS,
+                            scalarMode: ScalarMode.USE_POINT_FIELD_DATA,
+                        });
+                        const actor = vtkActor.newInstance();
+                        actor.getProperty().setColor(...type.rgb);
+                        actor.getProperty().setOpacity(opacity);
+                        actor.setMapper(mapper);
+                        if (shapeData.getClassName() == 'vtkPolyData') {
+                            mapper.setInputData(shapeData);
+                        } else if (cachedMarchingCubes.value[cacheLabel]) {
+                            mapper.setInputData(cachedMarchingCubes.value[cacheLabel])
+                        } else {
+                            const marchingCube = vtkImageMarchingCubes.newInstance({
+                                contourValue: 0.001,
+                                computeNormals: true,
+                                mergePoints: true,
+                            });
+                            marchingCube.setInputData(shapeData)
+                            mapper.setInputConnection(marchingCube.getOutputPort());
+                            cachedMarchingCubes.value[cacheLabel] = marchingCube.getOutputData()
+                        }
+                        if (showDifferenceFromMeanMode.value) {
+                            this.showDifferenceFromMean(mapper, renderer, label, domainIndex)
+                        }
+                        renderer.addActor(actor);
+                    }
+                )
             }
         )
     },
-    showDifferenceFromMean(mapper, renderer) {
-        if (!analysisFileShown.value
-            || !currentAnalysisFileParticles.value
-            || !meanAnalysisFileParticles.value
-            || !this.metaData.current
-            || !this.metaData.mean) return
+    showDifferenceFromMean(mapper, renderer, label, index) {
+        if (!analysisFilesShown.value
+            || !currentAnalysisParticlesFiles.value
+            || !meanAnalysisParticlesFiles.value
+            || !this.metaData[label]
+            || !this.metaData[label][index]) return
+
         // color values should be between 0 and 1
         // 0.5 is green, representing no difference between particles
-        const currentPoints = this.metaData.current.points.getPoints().getData()
-        const meanPoints = this.metaData.mean.points.getPoints().getData()
+        const targetMetaData = this.metaData[label][index]
+        const currentPoints = targetMetaData.current.points.getPoints().getData()
+        const meanPoints = targetMetaData.mean.points.getPoints().getData()
 
-        const particleComparisonKey = `${currentAnalysisFileParticles.value}_${meanAnalysisFileParticles.value}`
+        const particleComparisonKey = currentAnalysisParticlesFiles.value[index]
         let colorValues;
         let vectorValues;
         if (
@@ -382,14 +386,14 @@ export default {
         );
 
         this.addShapes(renderer, label, shapes.map(({ shape }) => shape));
-        shapes.map(({ points }) => points).forEach((pointSet) => {
+        shapes.map(({ points }) => points).forEach((pointSet, i) => {
             if (pointSet.getNumberOfPoints() > 0) {
-                this.addPoints(renderer, pointSet);
+                this.addPoints(renderer, pointSet, i);
             }
         })
-        shapes.map(({ landmarks }) => landmarks).forEach((landmarkSet) => {
+        shapes.map(({ landmarks }) => landmarks).forEach((landmarkSet, i) => {
             if (landmarkSet && landmarkSet.getNumberOfPoints() > 0) {
-                this.addPoints(renderer, landmarkSet, true);
+                this.addPoints(renderer, landmarkSet, i, true);
             }
         })
 
diff --git a/web/shapeworks/src/store/index.ts b/web/shapeworks/src/store/index.ts
index 183b53de..a170da47 100644
--- a/web/shapeworks/src/store/index.ts
+++ b/web/shapeworks/src/store/index.ts
@@ -45,11 +45,11 @@ export const groomedShapesForOriginalDataObjects = ref<Record<string, Record<num
 
 export const analysis = ref<Analysis>();
 
-export const analysisFileShown = ref<string>();
+export const analysisFilesShown = ref<string[]>();
 
-export const currentAnalysisFileParticles = ref<string>();
+export const currentAnalysisParticlesFiles = ref<string[]>();
 
-export const meanAnalysisFileParticles = ref<string>();
+export const meanAnalysisParticlesFiles = ref<string[]>();
 
 export const layersShown = ref<string[]>(["Original"])
 
diff --git a/web/shapeworks/src/store/methods.ts b/web/shapeworks/src/store/methods.ts
index 69061952..24dc842d 100644
--- a/web/shapeworks/src/store/methods.ts
+++ b/web/shapeworks/src/store/methods.ts
@@ -15,10 +15,20 @@ import {
      reconstructionsForOriginalDataObjects,
      cachedParticleComparisonColors,
      cachedParticleComparisonVectors,
-     meanAnalysisFileParticles,
      allDatasets,
      goodBadAngles,
-     landmarkColorList, landmarkInfo,
+     landmarkColorList,
+     landmarkInfo,
+     analysisExpandedTab,
+     selectedDataObjects,
+     particleSize,
+     analysisFilesShown,
+     meanAnalysisParticlesFiles,
+     currentAnalysisParticlesFiles,
+     goodBadMaxAngle,
+     showDifferenceFromMeanMode,
+     showGoodBadParticlesMode,
+     analysisAnimate,
 } from ".";
 import imageReader from "@/reader/image";
 import pointsReader from "@/reader/points";
@@ -39,6 +49,32 @@ import { layers, COLORS } from "./constants";
 import { getDistance, hexToRgb } from "@/helper";
 import router from "@/router";
 
+export const resetState = () => {
+    selectedDataObjects.value = [];
+    layersShown.value = ['Original'];
+    landmarkInfo.value = undefined;
+    landmarkColorList.value = [];
+    currentTasks.value = {};
+    jobProgressPoll.value = undefined;
+    particleSize.value = 2;
+    analysis.value = undefined;
+    analysisExpandedTab.value = 0;
+    analysisFilesShown.value = undefined;
+    currentAnalysisParticlesFiles.value = undefined;
+    meanAnalysisParticlesFiles.value = undefined;
+    goodBadAngles.value = undefined;
+    goodBadMaxAngle.value = 45;
+    showGoodBadParticlesMode.value = false;
+    showDifferenceFromMeanMode.value = false;
+    cachedMarchingCubes.value = {};
+    cachedParticleComparisonVectors.value = {};
+    cachedParticleComparisonColors.value = {};
+    landmarkInfo.value = undefined;
+    landmarkColorList.value = [];
+    analysisExpandedTab.value = 0;
+    analysisAnimate.value = false;
+}
+
 export const loadDataset = async (datasetId: number) => {
     // Only reload if something has changed
     if (selectedDataset.value?.id != datasetId) {
@@ -79,10 +115,10 @@ export const loadProjectsForDataset = async (datasetId: number) => {
 
 export const selectProject = (projectId: number | undefined) => {
     if (projectId) {
+        resetState();
         selectedProject.value = allProjectsForDataset.value.find(
             (project: Project) => project.id == projectId,
         )
-        layersShown.value = ["Original"]
         getLandmarks();
     }
 }
@@ -141,10 +177,19 @@ export async function spawnJob(action: string, payload: Record<string, any>): Pr
     if (!projectId) return undefined
     switch (action) {
         case 'groom':
+            layersShown.value = layersShown.value.filter(
+                (l) => l !== 'Groomed'
+            )
             return (await groomProject(projectId, payload))?.data
         case 'optimize':
+            layersShown.value = layersShown.value.filter(
+                (l) => l !== 'Particles'
+            )
             return (await optimizeProject(projectId, payload))?.data
         case 'analyze':
+            layersShown.value = layersShown.value.filter(
+                (l) => l !== 'Reconstructed'
+            )
             return (await analyzeProject(projectId, payload as AnalysisParams))?.data
         default:
             break;
@@ -224,7 +269,10 @@ export async function fetchJobResults(taskName: string) {
                 ([cachedLabel]) => !cachedLabel.includes(layerName)
             )
         )
-        if (!layersShown.value.includes(layerName)) layersShown.value.push(layerName)
+        const layer = layers.value.find((l) => l.name === layerName)
+        if (layer?.available() && !layersShown.value.includes(layerName)) {
+            layersShown.value = [...layersShown.value, layerName]
+        }
     }
 }
 
@@ -299,36 +347,46 @@ export function cacheComparison(colorValues: number[], vectorValues: number[][],
     cachedParticleComparisonVectors.value[particleComparisonKey] = vectorValues
 }
 
-export async function cacheAllComparisons(comparisons: CacheComparison[]) {
+export async function cacheAllComparisons(comparisons: CacheComparison[][]) {
     if (comparisons !== undefined) {
         const cachePrep = await Promise.all(comparisons?.map(async (g) => {
-            const particleComparisonKey = `${g.particles}_${meanAnalysisFileParticles.value}`;
+            return await Promise.all(
+                g.map(async (domain, i) => {
+                    if(meanAnalysisParticlesFiles.value && meanAnalysisParticlesFiles.value.length > i){
+                        const particleComparisonKey = domain.particles;
 
-            if (!cachedParticleComparisonColors.value[particleComparisonKey]) { // if the comparison is NOT already cached
-                const compareToPoints = await pointsReader(meanAnalysisFileParticles.value);
-                const currentPoints = await pointsReader(g.particles);
+                        if (!cachedParticleComparisonColors.value[particleComparisonKey]) { // if the comparison is NOT already cached
+                            const compareToPoints = await pointsReader(meanAnalysisParticlesFiles.value[i]);
+                            const currentPoints = await pointsReader(domain.particles);
 
-                const currentMesh = await imageReader(g.file, "current_mesh.vtk");
+                            const currentMesh = await imageReader(domain.file, "current_mesh.vtk");
 
-                return {
-                    "compareTo": {
-                        points: compareToPoints.getPoints().getData(),
-                        particleUrl: meanAnalysisFileParticles.value,
-                    },
-                    "current":  {
-                        points: currentPoints.getPoints().getData(),
-                        mapper: generateMapper(currentMesh),
-                        particleUrl: g.particles,
-                    },
-                }
-            }
+                            return {
+                                "compareTo": {
+                                    points: compareToPoints.getPoints().getData(),
+                                    particleUrl: domain,
+                                },
+                                "current":  {
+                                    points: currentPoints.getPoints().getData(),
+                                    mapper: generateMapper(currentMesh),
+                                    particleUrl: domain.particles,
+                                },
+                            }
+                        }
+                    }
+                })
+            )
         }))
 
         cachePrep.forEach((g) => {
             if (g !== undefined) {
-                const { current, compareTo } = g;
-                const comparisons = calculateComparisons(current.mapper, current.points as number[], compareTo.points as number[])
-                cacheComparison(comparisons.colorValues, comparisons.vectorValues, `${current.particleUrl}_${compareTo.particleUrl}`);
+                g.forEach((c) => {
+                    if (c){
+                        const { current, compareTo } = c;
+                        const comparisons = calculateComparisons(current.mapper, current.points as number[], compareTo.points as number[])
+                        cacheComparison(comparisons.colorValues, comparisons.vectorValues, current.particleUrl);
+                    }
+                })
             }
         })
     }
diff --git a/web/shapeworks/src/types/index.ts b/web/shapeworks/src/types/index.ts
index 33dc87f4..1f2a2d32 100644
--- a/web/shapeworks/src/types/index.ts
+++ b/web/shapeworks/src/types/index.ts
@@ -64,6 +64,11 @@ export interface AnalysisGroup {
     ratio: number,
 }
 
+export interface AnalysisMeanShape {
+    file: string,
+    particles: string,
+}
+
 export interface Analysis {
     id: number,
     created: string,
@@ -71,8 +76,7 @@ export interface Analysis {
     charts: AnalysisChart[],
     good_bad_angles: number[][],
     groups: AnalysisGroup[],
-    mean_shape: string,
-    mean_particles: string;
+    mean_shapes: AnalysisMeanShape[],
     modes: AnalysisMode[],
 }
 
diff --git a/web/shapeworks/src/views/Main.vue b/web/shapeworks/src/views/Main.vue
index 5c950cc1..d93fc0a7 100644
--- a/web/shapeworks/src/views/Main.vue
+++ b/web/shapeworks/src/views/Main.vue
@@ -26,9 +26,9 @@ import {
     selectedProject,
     loadProjectForDataset,
     reconstructionsForOriginalDataObjects,
-    analysisFileShown,
-    meanAnalysisFileParticles,
-    currentAnalysisFileParticles,
+    analysisFilesShown,
+    meanAnalysisParticlesFiles,
+    currentAnalysisParticlesFiles,
     switchTab,
     landmarkColorList,
     jobAlreadyDone,
@@ -79,7 +79,7 @@ export default {
         const rows = ref<number>(1);
         const cols = ref<number>(1);
         const renderData = ref<Record<string, ShapeData[]>>({});
-        const renderMetaData = ref<Record<string, ShapeData>>({});
+        const renderMetaData = ref<Record<string, ShapeData[]>>({});
 
         onMounted(async () => {
             try {
@@ -104,7 +104,7 @@ export default {
 
         async function toSelectPage() {
             selectedProject.value = undefined;
-            analysisFileShown.value = undefined;
+            analysisFilesShown.value = undefined;
             if (allProjectsForDataset.value.length === 0 && selectedDataset.value) {
                 loadProjectsForDataset(selectedDataset.value.id);
             }
@@ -151,49 +151,51 @@ export default {
 
         async function refreshRender() {
             let newRenderData = {}
-            renderMetaData.value = {}
+            let newRenderMetaData = {}
             const groupedSelections: Record<string, DataObject[]> = groupBy(selectedDataObjects.value, 'subject')
-            if (analysisFileShown.value) {
-                const currParticles = await pointsReader(
-                    currentAnalysisFileParticles.value
-                )
-                const meanParticles = await pointsReader(
-                    meanAnalysisFileParticles.value
-                )
-                renderMetaData.value = {
-                    "mean": {
-                        shape: await imageReader(undefined),
-                        points: meanParticles,
-                    },
-                    "current": {
-                        shape: await imageReader(undefined),
-                        points: currParticles,
-                    }
-                }
-                newRenderData =
-                    (analysisExpandedTab.value === 0) ?
-                    {
-                        "PCA": [{
-                            shape: await imageReader(
-                                analysisFileShown.value,
-                                shortFileName(analysisFileShown.value),
-                            ),
-                            points: await pointsReader(
-                                currentAnalysisFileParticles.value
+
+            if (
+                analysisFilesShown.value?.length &&
+                currentAnalysisParticlesFiles.value?.length &&
+                meanAnalysisParticlesFiles.value?.length
+            ) {
+                // populate newRenderData and renderMetaData
+                await Promise.all(
+                    analysisFilesShown.value.map(async (f, i) => {
+                        if(meanAnalysisParticlesFiles.value && currentAnalysisParticlesFiles.value){
+                            const key = analysisExpandedTab.value === 0 ? "PCA": "GROUP"
+                            if (!newRenderData[key]) {
+                                newRenderData[key] = []
+                            }
+                            newRenderData[key].push({
+                                shape: [  // only one layer, no overlapping shapes
+                                    await imageReader(f, shortFileName(f))
+                                ],
+                                points: await pointsReader(currentAnalysisParticlesFiles.value[i])
+                            })
+
+                            if (!newRenderMetaData[key]) {
+                                newRenderMetaData[key] = []
+                            }
+                            const meanParticles = await pointsReader(
+                                meanAnalysisParticlesFiles.value[i]
                             )
-                        }]
-                    }:
-                    {
-                        "GROUP": [{
-                            shape: await imageReader(
-                                analysisFileShown.value,
-                                shortFileName(analysisFileShown.value),
-                            ),
-                            points: await pointsReader(
-                                currentAnalysisFileParticles.value
+                            const currParticles = await pointsReader(
+                                currentAnalysisParticlesFiles.value[i]
                             )
-                        }]
-                    }
+                            newRenderMetaData[key].push({
+                                "mean": {
+                                    shape: await imageReader(undefined),
+                                    points: meanParticles,
+                                },
+                                "current": {
+                                    shape: await imageReader(undefined),
+                                    points: currParticles,
+                                }
+                            })
+                        }
+                    })
+                )
             } else {
                 newRenderData = Object.fromEntries(
                     await Promise.all(Object.entries(groupedSelections).map(
@@ -297,6 +299,7 @@ export default {
                 rows.value = numGroups;
             }
             renderData.value = newRenderData
+            renderMetaData.value = newRenderMetaData
         }
 
         const debouncedRefreshRender = _.debounce(refreshRender, 300)
@@ -304,9 +307,9 @@ export default {
         watch(drawer, prepareDrawer)
         watch(selectedDataObjects, debouncedRefreshRender)
         watch(layersShown, debouncedRefreshRender)
-        watch(analysisFileShown, debouncedRefreshRender)
+        watch(analysisFilesShown, debouncedRefreshRender, {deep: true})
         watch(landmarkColorList, debouncedRefreshRender)
-        watch(meanAnalysisFileParticles, debouncedRefreshRender)
+        watch(meanAnalysisParticlesFiles, debouncedRefreshRender, {deep: true})
         watch(tab, switchTab)
 
         return {
@@ -326,7 +329,7 @@ export default {
             refreshRender,
             jobAlreadyDone,
             particleSize,
-            analysisFileShown,
+            analysisFilesShown,
         }
     }
 }
@@ -411,7 +414,7 @@ export default {
         </v-card>
 
         <v-card :style="renderAreaStyle" class="render-area pa-3">
-            <template v-if="selectedDataObjects.length > 0 || analysisFileShown">
+            <template v-if="selectedDataObjects.length > 0 || analysisFilesShown && analysisFilesShown.length">
                 <shape-viewer
                     :data="renderData"
                     :metaData="renderMetaData"