Haha, it builds! God only knows if it works though.

src/Segmentation/SeededRegionGrower.cxx

@@ -1,47 +1,87 @@
 #include "SeededRegionGrower.h"
 
 #include <algorithm>
+#include <climits>
 #include <vector>
 #include <set>
 #include <map>
 
 class RatedVox {
-  int delta;
-  DICOMImage::IndexType idx;
+  public:
+    int delta;
+    DICOMImage::IndexType idx;
 
-  RatedVox(DICOMImage::IndexType index) {
-    this->idx = index;
-  }
+    RatedVox(DICOMImage::IndexType index) {
+      this->idx = index;
+    }
+
+    RatedVox(DICOMImage::IndexType index, DICOMImageP image, Region *region) {
+      this->idx = index;
+      this->delta = RatedVox::ComputeDelta(index, image, region);
+    }
+
+    IndexList GetNeighbors(DICOMImageP image) {
+      return SeededRegionGrower::GetNeighbors(image, this->idx);
+    }
+
+    static short ComputeDelta(
+        DICOMImage::IndexType idx, DICOMImageP image, Region *region) { 
+      PixelType intensity = image->GetPixel(idx);
+
+      return std::abs(intensity - region->mean);
+    }
+};
+
+struct VoxComp {
+ bool operator()(const RatedVox *v1, const RatedVox *v2) const {
+   return v1->delta < v2->delta;
+ }
+};
+
+bool IndexComp::operator()(
+   const DICOMImage::IndexType& i1, const DICOMImage::IndexType& i2) const {
+ for (int i=0; i<3; i++) {
+   if (i1[i] < i2[i]) {
+     return true;
+   } else if (i1[i] > i2[i]) {
+     return false;
+   }
+ }
+
+ // If i1 and i2 are equivalent
+ return false;
 }
 
-class Region {
-  string name;
-  PixelType mean;
-  IndexSet members;
+Region::Region(std::string name, IndexList seeds, DICOMImageP image) {
+  this->name = name;
 
-  Region(string name, RegionMembers seeds, DICOMImageP image) {
-    int sum = 0,
-      i = 0;
+  for (IndexList::iterator it=seeds.begin(); it!=seeds.end(); it++) {
+    DICOMImage::IndexType idx = *it;
+    this->AddPixel(idx, image);
+  }
+};
 
-    this->name = name;
+void Region::AddPixel(DICOMImage::IndexType idx, DICOMImageP image) {
+  this->members.insert(idx);
 
-    for(RegionMembers::iterator it=region_seeds.begin();
-        it != it.end(); it++) {
-      DICOMImage::IndexType idx = it*;
-      sum += image->GetPixel(idx);
-      i++;
+  this->sum += image->GetPixel(idx);
+  this->count++;
+  this->mean = this->sum / this->count;
+}
 
-      this->IndexSet.insert(idx);
-    }
 
-    this->mean = sum / i;
-  };
+bool Region::IsMember(DICOMImage::IndexType idx) {
+  return !!this->members.count(idx);
 }
 
-struct VoxComp {
- bool operator()(const RatedVox& v1, const RatedVox& v2) {
-   return v1.delta < v2.delta && v1.delta == v2.delta;
- }
+void SeededRegionGrower::FilterTouched(
+    IndexSet &touched, IndexList *idx_list) {
+  for (IndexList::iterator it=idx_list->begin(); it!=idx_list->end(); it++) {
+    if (touched.count(*it)) {
+      it = idx_list->erase(it);
+      it--;
+    }
+  }
 }
 
 DICOMImage::IndexType SeededRegionGrower::ConvertOffset(
@@ -53,9 +93,9 @@ DICOMImage::IndexType SeededRegionGrower::ConvertOffset(
   return idx;
 }
 
-RegionMembers SeededRegionGrower::GetNeighbors(
+IndexList SeededRegionGrower::GetNeighbors(
     DICOMImageP image, DICOMImage::IndexType idx) {
-  RegionMembers neighbors;
+  IndexList neighbors;
   DICOMImage::RegionType bounds = image->GetLargestPossibleRegion();
 
   for (int x=-1; x<2; x++) {
@@ -69,7 +109,7 @@ RegionMembers SeededRegionGrower::GetNeighbors(
                                            idx[1] + y,
                                            idx[2] + z }};
 
-        if (bounds->IsInside(new_idx)) {
+        if (bounds.IsInside(new_idx)) {
           neighbors.push_front(new_idx);
         }
       }
@@ -79,6 +119,31 @@ RegionMembers SeededRegionGrower::GetNeighbors(
   return neighbors;
 }
 
+Region *SeededRegionGrower::GetBestBorderingRegion(
+    std::set<Region*> &regions, DICOMImage::IndexType idx, DICOMImageP image) {
+  IndexList neighbors = SeededRegionGrower::GetNeighbors(image, idx);
+  Region *best_region;
+  PixelType best_delta = SHRT_MAX,
+    current_delta;
+
+  for (IndexList::iterator it=neighbors.begin(); it!=neighbors.end(); it++) {
+    for (std::set<Region*>::iterator it2=regions.begin();
+        it2!=regions.end(); ++it2) {
+      if ((*it2)->IsMember(*it)) {
+        current_delta = RatedVox::ComputeDelta(idx, image, *it2);
+
+        if (current_delta < best_delta) {
+          best_delta = current_delta;
+          best_region = *it2;
+        }
+      }
+    }
+  }
+
+  return best_region;
+}
+
+
 DICOMImage::Pointer SeededRegionGrower::LoadImage(std::string path) {
   ImageIOType::Pointer gdcmIO = ImageIOType::New();
 
@@ -98,30 +163,68 @@ DICOMImage::Pointer SeededRegionGrower::LoadImage(std::string path) {
 }
 
 DICOMImage::Pointer SeededRegionGrower::Segment(
-    DICOMImageP image, Regions seeds) {
-  Regions regions;
+    DICOMImageP image, RegionSeeds seeds) {
+  std::set<Region*> regions;
   IndexSet touched;
 
   // Not a list and only weakly sorted, w/e easier to read next to the '94 paper
-  std::vector<RatedVox> SSL;
+  std::vector<RatedVox*> SSL;
 
-  for (Regions::iterator it=seeds.begin(); it != seeds.end(); ++it) {
+  for (RegionSeeds::iterator it=seeds.begin(); it != seeds.end(); ++it) {
     // Copy seeds into regions (because each seed of a region is definitely
     // going to belong to it) and we'll grow as we iterate;
-    string region_name = it->first;
-    RegionMembers region_seeds = *(it->second);
+    std::string region_name = it->first;
+    IndexList region_seeds = *(it->second);
 
-    regions.insert(
-      std::pair<std::string, RegionMembers>(region_name, region_ssds));
+    Region *region = new Region(region_name, region_seeds, image);
+    regions.insert(region);
 
     // Insert each seed point into the touched set so we don't visit it
-    // and add it's neighbors to the SSL
-    for(RegionMembers::iterator it2=region_seeds.begin();
-        it2 != it2.end(); it2++) {
-      DICOMImage::IndexType idx = *it;
+    // and add them to the SSL so their neighbors get considered next iter
+    for(IndexList::iterator it2=region_seeds.begin();
+        it2 != region_seeds.end(); it2++) {
+      DICOMImage::IndexType idx = *it2;
       touched.insert(idx);
 
-      RegionMembers neighbors = SeededRegionGrower::GetNeighbors(image, idx);
+      RatedVox *seed = new RatedVox(idx, image, region);
+      region->AddPixel(idx, image);
+
+      // And add each seed point's neighbors to the SSL
+      IndexList neighbors = seed->GetNeighbors(image);
+      for (IndexList::iterator it3=neighbors.begin();
+          it3!=neighbors.end(); it3++) {
+        RatedVox *neighbor = new RatedVox(*it3, image, region);
+        SSL.push_back(neighbor);
+        std::push_heap(SSL.begin(), SSL.end(), VoxComp());
+      }
+    }
+  }
+
+  while (!SSL.empty()) {
+    // Vectors are straight retarded, std is a pile of junk
+    std::pop_heap(SSL.begin(), SSL.end(), VoxComp());
+    RatedVox *vox = SSL.back(),
+      *candidate;
+    SSL.pop_back();
+
+    Region *best_region = SeededRegionGrower::GetBestBorderingRegion(
+        regions, vox->idx, image);
+
+    best_region->AddPixel(vox->idx, image);
+
+    IndexList neighbors = vox->GetNeighbors(image),
+      neighbors2;
+    SeededRegionGrower::FilterTouched(touched, &neighbors);
+
+    for (IndexList::iterator it=neighbors.begin(); it!=neighbors.end(); it++) {
+      best_region = SeededRegionGrower::GetBestBorderingRegion(
+        regions, *it, image);
+
+      candidate = new RatedVox(*it, image, best_region);
+      SSL.push_back(candidate);
+      std::push_heap(SSL.begin(), SSL.end(), VoxComp());
     }
   }
+
+  return NULL;
 }

src/Segmentation/SeededRegionGrower.h

@@ -14,31 +14,64 @@
 #include "itkImageSeriesReader.h"
 
 typedef short PixelType;
-typedef std::list<DICOMImage::IndexType> RegionMembers;
-typedef std::set<DICOMImage::IndexType> IndexSet;
 typedef itk::Image<PixelType, 3> DICOMImage;
+typedef std::list<DICOMImage::IndexType> IndexList;
 typedef DICOMImage::Pointer DICOMImageP;
 typedef itk::ImageSeriesReader<DICOMImage> ReaderType;
 typedef itk::GDCMImageIO ImageIOType;
 typedef itk::GDCMSeriesFileNames InputNamesGeneratorType;
 typedef itk::NumericSeriesFileNames OutputNamesGeneratorType;
-typedef std::map<std::string, *RegionMembers> Regions;
-typedef std::map<std::string, std::list<DICOMImage::OffsetType> >
-  SpacialRegions;
+typedef std::map<std::string, IndexList*> RegionSeeds;
+
+struct IndexComp {
+ bool operator()(
+     const DICOMImage::IndexType &i1,
+     const DICOMImage::IndexType &i2) const;
+};
+
+typedef std::set<DICOMImage::IndexType, IndexComp> IndexSet;
+
+class Region {
+  public:
+    std::string name;
+    PixelType mean;
+  private:
+    int count, sum;
+    IndexSet members;
+
+  public:
+    Region(
+        std::string name,
+        IndexList seeds,
+        DICOMImageP image);
+    void AddPixel(
+        DICOMImage::IndexType idx,
+        DICOMImageP image);
+    bool IsMember(
+        DICOMImage::IndexType idx);
+};
 
 class SeededRegionGrower {
   public:
-    static DICOMImageP LoadImage(std::string path);
-    static DICOMImageP SegmentSpacial(
-        DICOMImageP image, SpacialRegions seeds);
+    static DICOMImageP LoadImage(
+        std::string path);
     static DICOMImageP Segment(
-        DICOMImageP image, Regions seeds);
+        DICOMImageP image,
+        RegionSeeds seeds);
     static DICOMImage::IndexType ConvertOffset(
-        DICOMImage::OffsetType offset, DICOMImage::SpacingType spacing);
-
+        DICOMImage::OffsetType offset,
+        DICOMImage::SpacingType spacing);
+    static IndexList GetNeighbors(
+        DICOMImageP image,
+        DICOMImage::IndexType idx);
+    static void FilterTouched(
+        IndexSet &touched,
+        IndexList *idx_list);
   private:
-    static RegionMembers GetNeighbors(
-      DICOMImageP image, DICOMImage::IndexType idx);
+    static Region *GetBestBorderingRegion(
+        std::set<Region*> &regions,
+        DICOMImage::IndexType idx,
+        DICOMImageP image);
 };
 
 #endif