diff --git a/bitmap_image.hpp b/bitmap_image.hpp
index 13ac766..702ead7 100644
--- a/bitmap_image.hpp
+++ b/bitmap_image.hpp
@@ -318,18 +318,37 @@ class bitmap_image
       return true;
    }
 
-   void reflective_image(bitmap_image& image)
+   void reflective_image(bitmap_image& image, const bool include_diagnols = false)
    {
       image.setwidth_height(3 * width_, 3 * height_,true);
-      image.copy_from(*this,width_,height_);
+
+      image.copy_from(*this, width_, height_);
+
       vertical_flip();
-      image.copy_from(*this,width_,0);
-      image.copy_from(*this,width_,2 * height_);
+
+      image.copy_from(*this, width_, 0);
+      image.copy_from(*this, width_, 2 * height_);
+
       vertical_flip();
       horizontal_flip();
-      image.copy_from(*this,0,height_);
-      image.copy_from(*this,2 * width_,height_);
+
+      image.copy_from(*this, 0, height_);
+      image.copy_from(*this, 2 * width_, height_);
+
       horizontal_flip();
+
+      if (include_diagnols)
+      {
+         bitmap_image tile = *this;
+
+         tile.vertical_flip();
+         tile.horizontal_flip();
+
+         image.copy_from(tile, 0 , 0);
+         image.copy_from(tile, 2 * width_, 0);
+         image.copy_from(tile, 2 * width_, 2 * height_);
+         image.copy_from(tile, 0, 2 * height_);
+      }
    }
 
    inline unsigned int width() const
@@ -1529,28 +1548,28 @@ class bitmap_image
    std::vector<unsigned char> data_;
 };
 
-struct rgb_store
+struct rgb_t
 {
-   unsigned char red;
+   unsigned char   red;
    unsigned char green;
-   unsigned char blue;
+   unsigned char  blue;
 };
 
-bool operator==(const rgb_store& c0, const rgb_store& c1)
+bool operator==(const rgb_t& c0, const rgb_t& c1)
 {
    return (c0.red   == c1  .red) &&
           (c0.green == c1.green) &&
           (c0.blue  == c1 .blue);
 }
 
-bool operator!=(const rgb_store& c0, const rgb_store& c1)
+bool operator!=(const rgb_t& c0, const rgb_t& c1)
 {
    return (c0.red   != c1  .red) ||
           (c0.green != c1.green) ||
           (c0.blue  != c1 .blue);
 }
 
-std::size_t hamming_distance(const rgb_store& c0, const rgb_store& c1)
+std::size_t hamming_distance(const rgb_t& c0, const rgb_t& c1)
 {
    std::size_t result = 0;
 
@@ -1561,6 +1580,62 @@ std::size_t hamming_distance(const rgb_store& c0, const rgb_store& c1)
    return result;
 }
 
+inline rgb_t make_colour(const unsigned int& red, const unsigned int& green, const unsigned int& blue)
+{
+   rgb_t result;
+
+   result.red   = static_cast<unsigned char>(red  );
+   result.green = static_cast<unsigned char>(green);
+   result.blue  = static_cast<unsigned char>(blue );
+
+   return result;
+}
+
+template <typename OutputIterator>
+inline void generate_colours(const std::size_t& steps, const rgb_t c0, const rgb_t& c1, OutputIterator out)
+{
+   double dr = ((double)c1.red   -  (double)c0.red   ) / steps;
+   double dg = ((double)c1.green -  (double)c0.green ) / steps;
+   double db = ((double)c1.blue  -  (double)c0.blue  ) / steps;
+
+   for (std::size_t i = 0; i < steps; ++i)
+   {
+      rgb_t c;
+
+      c.red   = static_cast<unsigned char>(c0.red   + (i * dr));
+      c.green = static_cast<unsigned char>(c0.green + (i * dg));
+      c.blue  = static_cast<unsigned char>(c0.blue  + (i * db));
+
+      *(out++) = c;
+   }
+}
+
+template <typename ResponseImage, typename Palette>
+inline std::size_t convert_rsp_to_image(const ResponseImage& resp_image, const Palette& palette, bitmap_image& image)
+{
+   if (
+        (resp_image.width () > image.width ()) ||
+        (resp_image.height() > image.height())
+      )
+      return 0;
+
+   for (std::size_t y = 0; y < resp_image.height(); ++y)
+   {
+      for (std::size_t x = 0; x < resp_image.width(); ++x)
+      {
+         const double v = resp_image(x,y);
+
+         unsigned int index = static_cast<unsigned int>((v < 0) ? 0 : v > (palette.size()) ? (palette.size() - 1) : v);
+
+         rgb_t c = palette[index];
+
+         image.set_pixel(x,y,c.red,c.green,c.blue);
+      }
+   }
+
+   return (resp_image.width() * resp_image.height());
+}
+
 inline void rgb_to_ycbcr(const unsigned int& length, double* red, double* green, double* blue,
                                                      double* y,   double* cb,    double* cr)
 {
@@ -1810,7 +1885,7 @@ inline void plasma(bitmap_image& image,
                    const double& c1,    const double& c2,
                    const double& c3,    const double& c4,
                    const double& roughness = 3.0,
-                   const rgb_store colormap[] = 0)
+                   const rgb_t colormap[] = 0)
 {
    // Note: c1,c2,c3,c4 -> [0.0,1.0]
 
@@ -1836,7 +1911,7 @@ inline void plasma(bitmap_image& image,
    }
    else
    {
-      rgb_store color = colormap[static_cast<unsigned int>(1000.0 * ((c1 + c2 + c3 + c4) / 4.0)) % 1000];
+      rgb_t color = colormap[static_cast<unsigned int>(1000.0 * ((c1 + c2 + c3 + c4) / 4.0)) % 1000];
       image.set_pixel(static_cast<unsigned int>(x),static_cast<unsigned int>(y),color.red,color.green,color.blue);
    }
 }
@@ -1889,7 +1964,7 @@ inline void hierarchical_psnr_r(const double& x,     const double& y,
                                 const bitmap_image& image1,
                                       bitmap_image& image2,
                                 const double& threshold,
-                                const rgb_store colormap[])
+                                const rgb_t colormap[])
 {
    if ((width <= 4.0) || (height <= 4.0))
    {
@@ -1903,7 +1978,7 @@ inline void hierarchical_psnr_r(const double& x,     const double& y,
 
       if (psnr < threshold)
       {
-         rgb_store c = colormap[static_cast<unsigned int>(1000.0 * (1.0 - (psnr / threshold)))];
+         rgb_t c = colormap[static_cast<unsigned int>(1000.0 * (1.0 - (psnr / threshold)))];
 
          image2.set_region(
                             static_cast<unsigned int>(x),
@@ -1926,7 +2001,7 @@ inline void hierarchical_psnr_r(const double& x,     const double& y,
    }
 }
 
-inline void hierarchical_psnr(bitmap_image& image1,bitmap_image& image2, const double threshold, const rgb_store colormap[])
+inline void hierarchical_psnr(bitmap_image& image1,bitmap_image& image2, const double threshold, const rgb_t colormap[])
 {
    if (
         (image1.width()  != image2.width ()) ||
@@ -2517,7 +2592,7 @@ class cartesian_canvas
    image_drawer draw_;
 };
 
-rgb_store convert_wave_length_nm_to_rgb(const double wave_length_nm)
+rgb_t convert_wave_length_nm_to_rgb(const double wave_length_nm)
 {
    // Credits: Dan Bruton http://www.physics.sfasu.edu/astro/color.html
    double red   = 0.0;
@@ -2572,7 +2647,7 @@ rgb_store convert_wave_length_nm_to_rgb(const double wave_length_nm)
    else
       factor = 0.0;
 
-   rgb_store result;
+   rgb_t result;
 
    const double gamma         =   0.8;
    const double intensity_max = 255.0;
@@ -2598,20 +2673,20 @@ double weighted_distance(const unsigned char r0, const unsigned char g0, const u
    return std::sqrt((diff_r * diff_r) + (diff_g * diff_g) + (diff_b * diff_b));
 }
 
-double weighted_distance(const rgb_store c0, const rgb_store c1)
+double weighted_distance(const rgb_t c0, const rgb_t c1)
 {
    return weighted_distance(c0.red, c0.green, c0.blue,
                             c1.red, c1.green, c1.blue);
 }
 
 template <typename Iterator>
-rgb_store find_nearest_color(const rgb_store& c, const Iterator begin, const Iterator end)
+rgb_t find_nearest_color(const rgb_t& c, const Iterator begin, const Iterator end)
 {
    if (0 == std::distance(begin,end))
       return c;
 
    double min_d     = std::numeric_limits<double>::max();
-   rgb_store result = *begin;
+   rgb_t result = *begin;
 
    for (Iterator itr = begin; itr != end; ++itr)
    {
@@ -2634,18 +2709,18 @@ rgb_store find_nearest_color(const rgb_store& c, const Iterator begin, const Ite
 
 template <template <typename,typename> class Sequence,
           typename Allocator>
-rgb_store find_nearest_color(const rgb_store& c, const Sequence<rgb_store,Allocator>& seq)
+rgb_t find_nearest_color(const rgb_t& c, const Sequence<rgb_t,Allocator>& seq)
 {
    return find_nearest_color(c, seq.begin(),seq.end());
 }
 
 template <std::size_t N>
-rgb_store find_nearest_color(const rgb_store& c, const rgb_store (&colors)[N])
+rgb_t find_nearest_color(const rgb_t& c, const rgb_t (&colors)[N])
 {
    return find_nearest_color(c, colors, colors + N);
 }
 
-double find_nearest_wave_length(const rgb_store& c, const double increment = 0.001)
+double find_nearest_wave_length(const rgb_t& c, const double increment = 0.001)
 {
    const double max_wave_length = 800.0; //800nm
 
@@ -2654,7 +2729,7 @@ double find_nearest_wave_length(const rgb_store& c, const double increment = 0.0
 
    for (double i = 0.0; i < max_wave_length; i += increment)
    {
-      rgb_store curr_rgb = convert_wave_length_nm_to_rgb(i);
+      rgb_t curr_rgb = convert_wave_length_nm_to_rgb(i);
       double    curr_d   = weighted_distance(c,curr_rgb);
 
       if (c == curr_rgb)
@@ -2672,6 +2747,64 @@ double find_nearest_wave_length(const rgb_store& c, const double increment = 0.0
    return min_wave_length;
 }
 
+template <typename T>
+class response_image
+{
+public:
+
+   response_image(const std::size_t& width, const std::size_t& height, const T null = T(0))
+   : width_ (width ),
+     height_(height),
+      null_ (null  )
+   {
+      data_.resize(width_ * height_);
+   }
+
+   std::size_t width () const { return  width_; }
+   std::size_t height() const { return height_; }
+
+   void set_all(const T& t)
+   {
+      std::fill_n(data_.begin(), data_.size(), t);
+   }
+
+   const T& operator()(const std::size_t& x, const std::size_t& y) const
+   {
+      if (y >= height_) return null_;
+      if (x >= width_ ) return null_;
+      return data_[width_ * y + x];
+   }
+
+   T& operator()(const std::size_t& x, const std::size_t& y)
+   {
+      if (y >= height_) return null_;
+      if (x >= width_ ) return null_;
+      return data_[width_ * y + x];
+   }
+
+   bool valid(const std::size_t& x, const std::size_t& y)
+   {
+      return ((x < width_ ) || (y < height_));
+   }
+
+   void inc_all(const T& v)
+   {
+      for (std::size_t i = 0; i < data_.size(); ++i)
+      {
+         data_[i] += v;
+      }
+   }
+
+private:
+
+   std::size_t    width_;
+   std::size_t    height_;
+   std::size_t    offset_x_;
+   std::size_t    offset_y_;
+   std::vector<T> data_;
+   T              null_;
+};
+
 enum palette_name
 {
   e_red,           e_scarlet,      e_vermilion,        e_tangelo,         e_orange,
@@ -2686,7 +2819,7 @@ enum palette_name
   e_raspberry,     e_crimson,      e_amaranth,         e_white,           e_black
 };
 
-const rgb_store palette_colormap[] = {
+const rgb_t palette_colormap[] = {
    {255,   0,   0}, {255,  31,   0}, {255,  63,   0}, {255,  95,   0}, {255, 127,   0},
    {255, 159,   0}, {255, 191,   0}, {255, 223,   0}, {255, 255,   0}, {223, 255,   0},
    {191, 255,   0}, {159, 255,   0}, {127, 255,   0}, { 95, 255,   0}, { 63, 255,   0},
@@ -2699,7 +2832,7 @@ const rgb_store palette_colormap[] = {
    {255,   0,  95}, {255,   0,  63}, {255,   0,  31}, {255, 255, 255}, {  0,   0,   0}
 };
 
-const rgb_store autumn_colormap[1000] = {
+const rgb_t autumn_colormap[1000] = {
    {255,   0,   0}, {255,   0,   0}, {255,   1,   0}, {255,   1,   0}, {255,   1,   0},
    {255,   1,   0}, {255,   2,   0}, {255,   2,   0}, {255,   2,   0}, {255,   2,   0},
    {255,   3,   0}, {255,   3,   0}, {255,   3,   0}, {255,   3,   0}, {255,   4,   0},
@@ -2902,7 +3035,7 @@ const rgb_store autumn_colormap[1000] = {
    {255, 254,   0}, {255, 254,   0}, {255, 254,   0}, {255, 255,   0}, {255, 255,   0}
 };
 
-const rgb_store copper_colormap[1000] = {
+const rgb_t copper_colormap[1000] = {
    {  0,   0,   0}, {  0,   0,   0}, {  1,   0,   0}, {  1,   1,   0}, {  1,   1,   1},
    {  2,   1,   1}, {  2,   1,   1}, {  2,   1,   1}, {  3,   2,   1}, {  3,   2,   1},
    {  3,   2,   1}, {  4,   2,   1}, {  4,   2,   2}, {  4,   3,   2}, {  4,   3,   2},
@@ -3105,7 +3238,7 @@ const rgb_store copper_colormap[1000] = {
    {255, 198, 126}, {255, 199, 126}, {255, 199, 127}, {255, 199, 127}, {255, 199, 127}
 };
 
-const rgb_store gray_colormap[1000] = {
+const rgb_t gray_colormap[1000] = {
    {255, 255, 255}, {255, 255, 255}, {254, 254, 254}, {254, 254, 254}, {254, 254, 254},
    {254, 254, 254}, {253, 253, 253}, {253, 253, 253}, {253, 253, 253}, {253, 253, 253},
    {252, 252, 252}, {252, 252, 252}, {252, 252, 252}, {252, 252, 252}, {251, 251, 251},
@@ -3308,7 +3441,7 @@ const rgb_store gray_colormap[1000] = {
    {  1,   1,   1}, {  1,   1,   1}, {  1,   1,   1}, {  0,   0,   0}, {  0,   0,   0}
 };
 
-const rgb_store hot_colormap[1000] = {
+const rgb_t hot_colormap[1000] = {
    { 11,   0,   0}, { 11,   0,   0}, { 12,   0,   0}, { 13,   0,   0}, { 13,   0,   0},
    { 14,   0,   0}, { 15,   0,   0}, { 15,   0,   0}, { 16,   0,   0}, { 17,   0,   0},
    { 17,   0,   0}, { 18,   0,   0}, { 19,   0,   0}, { 19,   0,   0}, { 20,   0,   0},
@@ -3511,7 +3644,7 @@ const rgb_store hot_colormap[1000] = {
    {255, 255, 251}, {255, 255, 252}, {255, 255, 253}, {255, 255, 254}, {255, 255, 255}
 };
 
-const rgb_store hsv_colormap[1000] = {
+const rgb_t hsv_colormap[1000] = {
    {255,   0,   0}, {255,   2,   0}, {255,   3,   0}, {255,   5,   0}, {255,   6,   0},
    {255,   8,   0}, {255,   9,   0}, {255,  11,   0}, {255,  12,   0}, {255,  14,   0},
    {255,  15,   0}, {255,  17,   0}, {255,  18,   0}, {255,  20,   0}, {255,  21,   0},
@@ -3714,7 +3847,7 @@ const rgb_store hsv_colormap[1000] = {
    {255,   0,  30}, {255,   0,  28}, {255,   0,  27}, {255,   0,  25}, {255,   0,  24}
 };
 
-const rgb_store jet_colormap[1000] = {
+const rgb_t jet_colormap[1000] = {
    { 29,   0, 102}, { 23,   0, 107}, { 17,   0, 112}, { 12,   0, 117}, {  6,   0, 122},
    {  0,   0, 127}, {  0,   0, 128}, {  0,   0, 129}, {  0,   0, 129}, {  0,   0, 130},
    {  0,   0, 131}, {  0,   0, 132}, {  0,   0, 133}, {  0,   0, 133}, {  0,   0, 134},
@@ -3917,7 +4050,7 @@ const rgb_store jet_colormap[1000] = {
    {122,   0,   9}, {117,   0,  18}, {112,   0,  27}, {107,   0,  36}, {102,   0,  45}
 };
 
-const rgb_store prism_colormap[1000] = {
+const rgb_t prism_colormap[1000] = {
    {255,   0,   0}, {255,   2,   0}, {255,   4,   0}, {255,   6,   0}, {255,   8,   0},
    {255,  10,   0}, {255,  11,   0}, {255,  13,   0}, {255,  15,   0}, {255,  17,   0},
    {255,  19,   0}, {255,  21,   0}, {255,  23,   0}, {255,  25,   0}, {255,  27,   0},
@@ -4120,7 +4253,7 @@ const rgb_store prism_colormap[1000] = {
    { 15, 255,   0}, { 11, 255,   0}, {  8, 255,   0}, {  4, 255,   0}, {  0, 255,   0}
 };
 
-const rgb_store vga_colormap[1000] = {
+const rgb_t vga_colormap[1000] = {
    {255, 255, 255}, {254, 254, 254}, {253, 253, 253}, {252, 252, 252}, {251, 251, 251},
    {250, 250, 250}, {249, 249, 249}, {248, 248, 248}, {247, 247, 247}, {246, 246, 246},
    {245, 245, 245}, {244, 244, 244}, {244, 244, 244}, {243, 243, 243}, {242, 242, 242},
@@ -4323,7 +4456,7 @@ const rgb_store vga_colormap[1000] = {
    {120,   0, 128}, {122,   0, 128}, {124,   0, 128}, {126,   0, 128}, {128,   0, 128}
 };
 
-const rgb_store yarg_colormap[1000] = {
+const rgb_t yarg_colormap[1000] = {
    {  0,   0,   0}, {  0,   0,   0}, {  1,   1,   1}, {  1,   1,   1}, {  1,   1,   1},
    {  1,   1,   1}, {  2,   2,   2}, {  2,   2,   2}, {  2,   2,   2}, {  2,   2,   2},
    {  3,   3,   3}, {  3,   3,   3}, {  3,   3,   3}, {  3,   3,   3}, {  4,   4,   4},
diff --git a/bitmap_test.cpp b/bitmap_test.cpp
index 21450eb..7b1ff8a 100644
--- a/bitmap_test.cpp
+++ b/bitmap_test.cpp
@@ -240,7 +240,7 @@ void test09()
    {
       for (unsigned int y = 0; y < dim; ++y)
       {
-         rgb_store col = jet_colormap[(x + y) % dim];
+         rgb_t col = jet_colormap[(x + y) % dim];
          image.set_pixel(x,y,col.red,col.green,col.blue);
       }
    }
@@ -471,17 +471,17 @@ void test18()
    {
       bitmap_image image(1000,180);
       image_drawer draw(image);
-      const rgb_store* colormap[9] = {
-                                       autumn_colormap,
-                                       copper_colormap,
-                                         gray_colormap,
-                                          hot_colormap,
-                                          hsv_colormap,
-                                          jet_colormap,
-                                        prism_colormap,
-                                          vga_colormap,
-                                         yarg_colormap
-                                     };
+      const rgb_t* colormap[9] = {
+                                   autumn_colormap,
+                                   copper_colormap,
+                                     gray_colormap,
+                                      hot_colormap,
+                                      hsv_colormap,
+                                      jet_colormap,
+                                    prism_colormap,
+                                      vga_colormap,
+                                     yarg_colormap
+                                 };
 
       for (unsigned int i = 0; i < image.width(); ++i)
       {
@@ -499,7 +499,7 @@ void test18()
       bitmap_image image(1000,500);
       image_drawer draw(image);
 
-      std::size_t palette_colormap_size = sizeof(palette_colormap) / sizeof(rgb_store);
+      std::size_t palette_colormap_size = sizeof(palette_colormap) / sizeof(rgb_t);
       std::size_t bar_width = image.width() / palette_colormap_size;
 
       for (std::size_t i = 0; i < palette_colormap_size; ++i)
@@ -536,7 +536,7 @@ void test19()
 
       for (std::size_t i = 0; i < 160; i++)
       {
-         std::size_t c_idx = i % (sizeof(palette_colormap) / sizeof(rgb_store));
+         std::size_t c_idx = i % (sizeof(palette_colormap) / sizeof(rgb_t));
 
          canvas.pen_color(palette_colormap[c_idx].red, palette_colormap[c_idx].green, palette_colormap[c_idx].blue);
 
@@ -565,7 +565,7 @@ void test19()
 
       for (double x = -500; x < 500; x += 3, ++i)
       {
-         std::size_t c_idx = i % (sizeof(palette_colormap) / sizeof(rgb_store));
+         std::size_t c_idx = i % (sizeof(palette_colormap) / sizeof(rgb_t));
 
          canvas.pen_color(palette_colormap[c_idx].red, palette_colormap[c_idx].green, palette_colormap[c_idx].blue);
 
@@ -582,12 +582,12 @@ void test19()
 
 void test20()
 {
-   const rgb_store* colormap[4] = {
-                                       hsv_colormap,
-                                       jet_colormap,
-                                     prism_colormap,
-                                       vga_colormap
-                                  };
+   const rgb_t* colormap[4] = {
+                                   hsv_colormap,
+                                   jet_colormap,
+                                 prism_colormap,
+                                   vga_colormap
+                              };
 
    const unsigned int fractal_width  = 1200;
    const unsigned int fractal_height =  800;
@@ -639,10 +639,10 @@ void test20()
                         (1000.0 * log2(1.75 + i - log2(log2(z))) / log2(max_iterations));
                      #undef log2
 
-                     rgb_store c0 = colormap[0][index];
-                     rgb_store c1 = colormap[1][index];
-                     rgb_store c2 = colormap[2][index];
-                     rgb_store c3 = colormap[3][index];
+                     rgb_t c0 = colormap[0][index];
+                     rgb_t c1 = colormap[1][index];
+                     rgb_t c2 = colormap[2][index];
+                     rgb_t c3 = colormap[3][index];
 
                      fractal_hsv  .set_pixel(x, y, c0.red, c0.green, c0.blue);
                      fractal_jet  .set_pixel(x, y, c1.red, c1.green, c1.blue);
@@ -701,10 +701,10 @@ void test20()
                   {
                      unsigned int index = static_cast<int>((1000.0 * i) / max_iterations);
 
-                     rgb_store c0 = colormap[0][index];
-                     rgb_store c1 = colormap[1][index];
-                     rgb_store c2 = colormap[2][index];
-                     rgb_store c3 = colormap[3][index];
+                     rgb_t c0 = colormap[0][index];
+                     rgb_t c1 = colormap[1][index];
+                     rgb_t c2 = colormap[2][index];
+                     rgb_t c3 = colormap[3][index];
 
                      fractal_hsv  .set_pixel(x, y, c0.red, c0.green, c0.blue);
                      fractal_jet  .set_pixel(x, y, c1.red, c1.green, c1.blue);
diff --git a/readme.md b/readme.md
index 525b3f5..7c7f2a8 100644
--- a/readme.md
+++ b/readme.md
@@ -1,5 +1,7 @@
 #### Description
-The C++ Bitmap Library consists of simple, robust, optimized and portable processing routines for the **24-bit per pixel** bitmap image format.
+The C++ Bitmap Library consists of simple, robust, optimized and
+portable processing routines for the **24-bit per pixel** bitmap
+image format.
 
 
 #### Capabilities
@@ -33,9 +35,9 @@ The C++ Bitmap Library implementation is compatible with the following C++ compi
 ----
 
 #### Simple Example 1
-The following example will open a bitmap image called *'input.bmp'* and count the
-number of pixels that have a red channel value of 111 and larger, then proceed
-to print the count to stdout.
+The following example will open a bitmap image called *'input.bmp'* and count
+the number of pixels that have a red channel value of 111 and larger, then
+proceed to print the count to stdout.
 
 ```c++
 #include <cstdio>
@@ -113,7 +115,8 @@ int main()
 ----
 
 #### Simple Example 3
-The following example will render the Mandelbrot set fractal and save the generated bitmap as *'mandelbrot_set.bmp'*.
+The following example will render the Mandelbrot set fractal and save
+the generated bitmap as *'mandelbrot_set.bmp'*.
 
 ```c++
 #include <cmath>
@@ -326,10 +329,11 @@ int main()
 ![ScreenShot](http://www.partow.net/programming/bitmap/images/lens_effect.png?raw=true "C++ Bitmap Library Magnifying Lens Effect Example - By Arash Partow")
 
 #### Simple Example 6
-The following example will render a baseline image using a combination of plasma and
-checkered pattern effects. Then proceed to apply a swirl distortion upon the base image.
-Finally both the base and the swirl distorted versions of the images will be saved to
-file as *'base.bmp'* and *'swirl_effect.bmp'* respectively.
+The following example will render a baseline image using a combination
+of plasma and checkered pattern effects. Then proceed to apply a swirl
+distortion upon the base image. Finally both the base and the swirl
+distorted versions of the images will be saved to file as *'base.bmp'*
+and *'swirl_effect.bmp'* respectively.
 
 ```c++
 #include <algorithm>
@@ -408,3 +412,165 @@ int main()
 ```
 
 ![ScreenShot](http://www.partow.net/programming/bitmap/images/swirl_effect.png?raw=true "C++ Bitmap Library Swirl Effect Example - By Arash Partow")
+
+#### Simple Example 7
+The following example will render a group of fireballs placed
+equidistant to their immediate neighbours following a Lissajous curve.
+The fireballs will then proceed move within the plane using the curve
+as their path. The example demonstrates the construction of a
+piecewise colour palette and the response_image functionality. After
+having 'simulated' N-frames, the final frame will be converted to a
+bitmap and then saved to file as *'fireballs.bmp'*.
+
+```c++
+#include <cmath>
+#include <iterator>
+#include <vector>
+#include "bitmap_image.hpp"
+
+struct lissajous_curve
+{
+   lissajous_curve(const double xs, const double ys)
+   : scale_x(xs),
+     scale_y(ys)
+   {}
+
+   inline double x(const double t) const { return scale_x * std::sin(    t); }
+   inline double y(const double t) const { return scale_y * std::cos(3 * t); }
+
+   double scale_x;
+   double scale_y;
+};
+
+int main()
+{
+   bitmap_image image(600,400);
+
+   image.clear();
+
+   lissajous_curve curve(image.width () / 2.0 - 25, image.height() / 2.0 - 25);
+
+   const double pi_   = 3.1415926535897932384626433832795028841971;
+   const double delta = (2.0 * pi_) / 100000.0;
+
+   const std::size_t max_fire_balls   = 30;
+   const std::size_t number_of_frames = 3000;
+
+   const double cooling_factor = 0.940; // [0,1]
+
+   // Arc-length of curve: x(t) = a0 * sin(4t), y(t) = a1 * cos(3t)
+   const double curve_length = 5102.0;
+
+   double segment_length  = curve_length / max_fire_balls;
+   double curr_seg_length = 0;
+
+   double prev_x = curve.x(0);
+   double prev_y = curve.y(0);
+
+   std::vector<double> fire_ball;
+
+   // Set the initial location for each fireball
+   for (double t = delta;  fire_ball.size() < max_fire_balls; t += delta)
+   {
+      double center_x = curve.x(t);
+      double center_y = curve.y(t);
+
+      double dx = (prev_x - center_x);
+      double dy = (prev_y - center_y);
+
+      curr_seg_length += std::sqrt((dx * dx) + (dy * dy));
+
+      prev_x = center_x;
+      prev_y = center_y;
+
+      if (curr_seg_length >= segment_length)
+      {
+         curr_seg_length = 0.0;
+         fire_ball.push_back(t);
+      }
+   }
+
+   response_image<double> resp_image(image.width(),image.height(), -1.0);
+   response_image<double> fb_misses (image.width(),image.height(), -1.0);
+
+   resp_image.set_all(-999.0);
+   fb_misses .set_all(   0.0);
+
+   std::vector<rgb_t> fire_palette;
+
+   // Baseline colours used fire palette
+   rgb_t black  = make_colour(0,    0,  0);
+   rgb_t red    = make_colour(255,  0,  0);
+   rgb_t yellow = make_colour(255,255,  0);
+   rgb_t white  = make_colour(255,255,255);
+
+   // Setup the fire palette:
+   // Black (Coolest - 0) --> Red --> Yellow --> White (Hottest - 999)
+   generate_colours(334, black,    red, std::back_inserter(fire_palette));
+   generate_colours(333,   red, yellow, std::back_inserter(fire_palette));
+   generate_colours(333,yellow,  white, std::back_inserter(fire_palette));
+
+   for (std::size_t k = 0; k < number_of_frames; ++k)
+   {
+      fb_misses.inc_all(1);
+
+      // Render fireballs on response image
+      for (std::size_t i = 0; i < fire_ball.size(); ++i)
+      {
+         double fb_x = curve.x(fire_ball[i]) + image.width () / 2.0;
+         double fb_y = curve.y(fire_ball[i]) + image.height() / 2.0;
+
+         // Draw circles with radii in the range [1,10]
+         for (double t = 0;  t < (2 * pi_); t += (((2.0 * pi_) / 360)))
+         {
+            for (double r = 1; r <= 10; ++r)
+            {
+               std::size_t rx = static_cast<std::size_t>(r * std::sin(t) + fb_x);
+               std::size_t ry = static_cast<std::size_t>(r * std::cos(t) + fb_y);
+
+               double heat_distortion = 50.0 * std::cos(delta * i) + 50; // per-frame in [0,100]
+
+               resp_image(rx,ry) = fire_palette.size() * 0.8 + heat_distortion + (::rand() % 100);
+               fb_misses (rx,ry) = 0;
+            }
+         }
+
+         // Move fireball to its next location
+         fire_ball[i] += delta;
+      }
+
+      // Apply cooling process to the entire plane
+      for (std::size_t y = 1; y < resp_image.height() - 1; ++y)
+      {
+         for (std::size_t x = 1; x < resp_image.width() - 1; ++x)
+         {
+            double avg = (
+                           resp_image(x - 1, y - 1) + resp_image(x    , y - 1) +
+                           resp_image(x + 1, y - 1) + resp_image(x - 1, y    ) +
+                           resp_image(x + 1, y    ) + resp_image(x    , y + 1) +
+                           resp_image(x - 1, y + 1) + resp_image(x + 1, y + 1)
+                         ) / (7.0 + cooling_factor);
+
+            // Only allow cooler averages to be applied
+            if (avg > resp_image(x, y))
+               continue;
+
+            // More rapidly cool points that haven't seen fireballs in the last N-frames
+            if (fb_misses(x,y) > 2000)
+               avg *= 0.90 + ((::rand() % 10) / 100.0);
+
+            // Clamp average in the range [0,999]
+            resp_image(x,y) = ((avg < 0.0) ? 0.0 : ((avg > 999.0) ? 999.0 : avg));
+         }
+      }
+   }
+
+   convert_rsp_to_image(resp_image,fire_palette,image);
+
+   image.save_image("fireballs.bmp");
+
+   return 0;
+}
+```
+
+![ScreenShot](http://www.partow.net/programming/bitmap/images/fireballs.png?raw=true "C++ Bitmap Library Fire Balls Example - By Arash Partow")