diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index 45688c72a4fb6c5a2983262f86bff107dbfffe8c..c943294a29f2961e39b3bb89febe1f2967cde112 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -2,7 +2,8 @@
 """
 @author: Laura C. Kühle, Soraya Terrab (sorayaterrab)
 
-TODO: Vectorize _get_cells() in Boxplot method
+TODO: Vectorize _get_cells() in Boxplot method -> Done
+TODO: Restructure Boxplot method
 TODO: Introduce lower/upper extreme outliers in Boxplot
     (each cell is also checked for neighboring domains if existing)
 TODO: Determine max_value for Theoretical only over highest degree
@@ -388,87 +389,81 @@ class Boxplot(WaveletDetector):
             List of indices for all detected troubled cells.
 
         """
-        # indexed_coeffs = [[multiwavelet_coeffs[0, i], i]
-        #                   for i in range(self._mesh.num_grid_cells)]
         coeffs = multiwavelet_coeffs[0]
-        # print(coeffs.shape)
 
         if self._mesh.num_grid_cells < self._fold_len:
             self._fold_len = self._mesh.num_grid_cells
 
-        num_folds = self._mesh.num_grid_cells//self._fold_len
-        troubled_cells = []
-        # troubled_cells_new = []
-
-        for fold in range(num_folds):
-            # indexed_fold = np.array(indexed_coeffs)[self._folds[fold]]
-            # sorted_fold_old = indexed_fold[indexed_fold[:, 0].argsort()]
-
-            sorted_fold = sorted(coeffs[self._folds[fold]])
-            # print(sorted_fold == sorted_fold_old[:, 0])
-
-            boundary_index = self._fold_len//4
-            balance_factor = self._fold_len/4.0 - boundary_index
-
-            first_quartile = (1-balance_factor) \
-                * sorted_fold[boundary_index-1] \
-                + balance_factor * sorted_fold[boundary_index]
-            third_quartile = (1-balance_factor) \
-                * sorted_fold[3*boundary_index-1]\
-                + balance_factor * sorted_fold[3*boundary_index]
-
-            lower_bound = first_quartile \
-                - self._whisker_len * (third_quartile-first_quartile)
-            upper_bound = third_quartile \
-                + self._whisker_len * (third_quartile-first_quartile)
-
-            # Adjust outer fences if flag is set
-            if self._adjust_outer_fences:
-                global_mean = np.mean(abs(coeffs))
-                lower_bound = min(-global_mean, lower_bound)
-                upper_bound = max(global_mean, upper_bound)
-
-            # # Check for lower extreme outliers and add respective cells
-            # for cell in sorted_fold:
-            #     if cell[0] < lower_bound:
-            #         troubled_cells.append(int(cell[1]))
-            #     else:
-            #         break
-            #
-            # # Check for upper extreme outliers and add respective cells
-            # for cell in sorted_fold[::-1][:]:
-            #     if cell[0] > upper_bound:
-            #         troubled_cells.append(int(cell[1]))
-            #     else:
-            #         break
-
-            # Check for extreme outlier and add respective cells
-            for cell in self._folds[fold]:
-                if (coeffs[cell] > upper_bound) \
-                        or (coeffs[cell] < lower_bound):
-                    troubled_cells.append(int(cell))
-
-            # print(upper_bound, lower_bound)
-            # print(sorted_fold_new)
-            # print(type(sorted_fold_new))
-            # print(sorted_fold_new > upper_bound)
-            # print(sorted_fold_new < lower_bound)
-            # test =
-            # print(type(test), test)
-            # print(list(test), list(test[0]))
-
-            # troubled_cells_new += list(np.flatnonzero(np.logical_or(
-            #     sorted_fold_new > upper_bound,
-            #     sorted_fold_new < lower_bound)).astype(int))
-            # print(troubled_cells_new)
-
-        # troubled_cells_new = sorted(troubled_cells_new)
-        # print(troubled_cells_new)
-        # print(troubled_cells)
-        # print(sorted(troubled_cells) == sorted(troubled_cells_new))
-        # print(type(troubled_cells_new[0]), type(troubled_cells[0]))
-
-        return sorted(troubled_cells)
+        boundary_index = self._fold_len//4
+        balance_factor = self._fold_len/4.0 - boundary_index
+
+        folds = coeffs[self._folds]
+        folds.sort()
+        first_quartiles = (1-balance_factor) \
+            * folds[:, boundary_index-1] \
+            + balance_factor * folds[:, boundary_index]
+        third_quartiles = (1-balance_factor) \
+            * folds[:, 3*boundary_index-1]\
+            + balance_factor * folds[:, 3*boundary_index]
+
+        lower_bounds = first_quartiles - self._whisker_len * (
+                third_quartiles-first_quartiles)
+        upper_bounds = third_quartiles + self._whisker_len * (
+                third_quartiles-first_quartiles)
+
+        if self._adjust_outer_fences:
+            global_mean = np.mean(abs(coeffs))
+            lower_bounds[lower_bounds > -global_mean] = -global_mean
+            upper_bounds[upper_bounds < global_mean] = global_mean
+
+        troubled_cells_new = np.flatnonzero(np.logical_or(
+            coeffs < np.repeat(lower_bounds, self._fold_len),
+            coeffs > np.repeat(upper_bounds, self._fold_len))).tolist()
+
+        # num_folds = self._mesh.num_grid_cells//self._fold_len
+        # troubled_cells = []
+        #
+        # lower_bound = np.zeros(num_folds)
+        # upper_bound = np.zeros(num_folds)
+        #
+        # for fold in range(num_folds):
+        #     sorted_fold = sorted(coeffs[self._folds[fold]])
+        #
+        #     first_quartile = (1-balance_factor) \
+        #         * sorted_fold[boundary_index-1] \
+        #         + balance_factor * sorted_fold[boundary_index]
+        #     third_quartile = (1-balance_factor) \
+        #         * sorted_fold[3*boundary_index-1]\
+        #         + balance_factor * sorted_fold[3*boundary_index]
+        #
+        #     lower_bound[fold] = first_quartile \
+        #         - self._whisker_len * (third_quartile-first_quartile)
+        #     upper_bound[fold] = third_quartile \
+        #         + self._whisker_len * (third_quartile-first_quartile)
+        #
+        #     # Adjust outer fences if flag is set
+        #     if self._adjust_outer_fences:
+        #         global_mean = np.mean(abs(coeffs))
+        #         lower_bound[fold] = min(-global_mean, lower_bound[fold])
+        #         upper_bound[fold] = max(global_mean, upper_bound[fold])
+        #
+        #     # Check for extreme outlier and add respective cells
+        #     for cell in self._folds[
+        #                 fold, self._num_overlapping_cells:
+        #                 -self._num_overlapping_cells]:
+        #         if (coeffs[cell] > upper_bound[fold]) \
+        #                 or (coeffs[cell] < lower_bound[fold]):
+        #             troubled_cells.append(int(cell))
+        #
+        # same = np.all(sorted(troubled_cells) == troubled_cells_new)
+        # if not same:
+        #     print(np.all(lower_bounds == lower_bound),
+        #           np.all(upper_bounds == upper_bound))
+        #     print(sorted(troubled_cells))
+        #     print(troubled_cells_new)
+
+        return troubled_cells_new
+        # return sorted(troubled_cells)
 
 
 class Theoretical(WaveletDetector):