diff --git a/czi-format/czi-parser/czi_file.cpp b/czi-format/czi-parser/czi_file.cpp
index b31da178aa5a64724388bcbf35efa09573905b0a..c3da3c0b2e7c990834efa8eb0819c38798d43f33 100644
--- a/czi-format/czi-parser/czi_file.cpp
+++ b/czi-format/czi-parser/czi_file.cpp
@@ -5,7 +5,7 @@ bool CziFile::is_master_file() const
     return ((header.filePart == 0) && (header.fileGuid == header.masterFileGuid));
 }
 
-void print_segment_header(const SegmentHeader &segmentHeader)
+void CziFile::print_segment_header(const SegmentHeader &segmentHeader) const
 {
     printf("---------SegmentHeader---------\n");
     printf("%-25s %15s\n", "SID", segmentHeader.sId.c_str());
@@ -14,7 +14,7 @@ void print_segment_header(const SegmentHeader &segmentHeader)
     printf("-------------------------------\n");
 }
 
-const char *pixel_type_str(const PixelType px)
+const char *CziFile::pixel_type_str(const PixelType px) const
 {
 
     switch (px)
@@ -49,7 +49,7 @@ const char *pixel_type_str(const PixelType px)
     return nullptr;
 }
 
-const char *pyramid_type_str(const PyramidType pt)
+const char *CziFile::pyramid_type_str(const PyramidType pt) const
 {
     switch (pt)
     {
@@ -69,7 +69,7 @@ const char *pyramid_type_str(const PyramidType pt)
     return nullptr;
 }
 
-const char *compression_type_str(const CompressionType ct)
+const char *CziFile::compression_type_str(const CompressionType ct) const
 {
     switch (ct)
     {
@@ -95,7 +95,7 @@ const char *compression_type_str(const CompressionType ct)
     return nullptr;
 }
 
-const char *dimension_type_str(const Dimension &d)
+const char *CziFile::dimension_type_str(const Dimension &d) const
 {
     switch (d)
     {
@@ -134,35 +134,46 @@ const char *dimension_type_str(const Dimension &d)
     return nullptr;
 }
 
-std::string dimension_char(const Dimension &d)
+const char *CziFile::dimension_char(const Dimension &d) const
 {
+    std::string result;
     switch (d)
     {
     case Dimension::X:
-        return "X";
+        result = "X";
+        break;
     case Dimension::Y:
-        return "Y";
+        result = "Y";
+        break;
     case Dimension::C:
-        return "C";
+        result = "C";
+        break;
     case Dimension::Z:
-        return "Z";
+        result = "Z";
+        break;
     case Dimension::T:
-        return "T";
+        result = "T";
+        break;
     case Dimension::R:
-        return "R";
+        result = "R";
+        break;
     case Dimension::S:
-        return "S";
+        result = "S";
+        break;
     case Dimension::I:
-        return "I";
+        result = "I";
+        break;
     case Dimension::B:
-        return "B";
+        result = "B";
+        break;
     case Dimension::M:
-        return "M";
+        result = "M";
+        break;
     case Dimension::H:
-        return "H";
+        result = "H";
+        break;
     case Dimension::V:
-        return "V";
-
+        result = "V";
         break;
     default:
     {
@@ -170,10 +181,10 @@ std::string dimension_char(const Dimension &d)
         break;
     }
     }
-    return "";
+    return result.c_str();
 }
 
-std::string dimension_stack_str(const std::vector<DimensionEntryDV1> &dims, bool includeSize)
+const char *CziFile::dimension_stack_str(const std::vector<DimensionEntryDV1> &dims, bool includeSize) const
 {
 
     std::string result;
@@ -182,7 +193,7 @@ std::string dimension_stack_str(const std::vector<DimensionEntryDV1> &dims, bool
         result.append(dimension_char(entry.dimension) + (includeSize ? ("(" + std::to_string(entry.size) + ")-") : ""));
     }
 
-    return includeSize ? result.substr(0, result.length() - 1) : result;
+    return (includeSize ? result.substr(0, result.length() - 1) : result).c_str();
 }
 
 void CziFile::report_verbose() const
@@ -199,7 +210,7 @@ void CziFile::report_verbose() const
                (int)entry,
                pixel_type_str(subBlockDirectory.entries[entry].pixelType),
                compression_type_str(subBlockDirectory.entries[entry].compression),
-               dimension_stack_str(subBlockDirectory.entries[entry].dimensions, true).c_str(),
+               dimension_stack_str(subBlockDirectory.entries[entry].dimensions, true),
                pyramid_type_str(subBlockDirectory.entries[entry].pyramidType));
 
         // printf("-------SubBlock directory entry %i-------\n", (int)entry);
diff --git a/czi-format/czi-parser/czi_file.h b/czi-format/czi-parser/czi_file.h
index c03c841e9b4aec9d1bc0113faee572b45ec74aa4..5a29aaa00de09d5af856a6633cb060a48a9e55a3 100644
--- a/czi-format/czi-parser/czi_file.h
+++ b/czi-format/czi-parser/czi_file.h
@@ -7,8 +7,18 @@
 //#include "../image/matrix.h"
 
 //TODO: Handle multi-file scenarios.
-struct CziFile
+class CziFile
 {
+private:
+  void print_segment_header(const SegmentHeader &segmentHeader) const;
+  const char *pixel_type_str(const PixelType px) const;
+  const char *pyramid_type_str(const PyramidType pt) const;
+  const char *compression_type_str(const CompressionType ct) const;
+  const char *dimension_type_str(const Dimension &d) const;
+  const char *dimension_char(const Dimension &d) const;
+  const char *dimension_stack_str(const std::vector<DimensionEntryDV1> &dims, bool includeSize) const;
+
+public:
   // Path to loaded czi file.
   std::string fileName;
   FileHeaderSegment header;
diff --git a/czi-format/czi-parser/image/matrix.h b/czi-format/czi-parser/image/matrix.h
index d9f848b791e26ec3d63332b9d3cfc3cd33861dd8..9944ad58c4d1deac49529f820de7981e0cbdb056 100644
--- a/czi-format/czi-parser/image/matrix.h
+++ b/czi-format/czi-parser/image/matrix.h
@@ -3,15 +3,14 @@
 
 #pragma once
 
-template <typename T>
 class Matrix
 {
   private:
     uint rowCount;
     uint colCount;
-    std::vector<T> data;
+    std::vector<BasePixel *> data;
 
-    bool equals(const Matrix<T> &other)
+    bool equals(const Matrix &other)
     {
         if ((this->rowCount != other.rows()) || (this->colCount != other.cols()))
             return false;
@@ -19,14 +18,20 @@ class Matrix
         size_t dataSize = this->data.size();
         for (size_t i = 0; i < dataSize; i++)
         {
-            //if (px_eq(data[i], other.data[i]))
-            if (this->data[i] != other.data[i])
+            if (data[i])
+            {
+                if (!this->data[i]->equals(other.data[i]))
+                    return false;
+            }
+            else if (other.data[i])
+            {
                 return false;
+            }
         }
         return true;
     }
 
-    std::function<T(std::vector<byte> &, ulong &)> get_pixel_parse_function(const PixelType pixelType)
+    std::function<BasePixel(std::vector<byte> &, ulong &)> get_pixel_parse_function(const PixelType pixelType)
     {
         switch (pixelType)
         {
@@ -36,8 +41,8 @@ class Matrix
         //     return "Gray16";
         // case PixelType::Gray32Float:
         //     return "Gray32Float";
-        case PixelType::Bgr24:
-            return bytes_to_bgr24pixel;
+        // case PixelType::Bgr24:
+        //     return bytes_to_bgr24pixel;
         // case PixelType::Bgr48:
         //     return "Bgr48";
         // case PixelType::Bgr96Float:
@@ -50,8 +55,8 @@ class Matrix
         //     return "Bgr192ComplexFloat";
         // case PixelType::Gray32:
         //     return "Gray32";
-        // case PixelType::Gray64:
-        //     return "Gray64";
+        case PixelType::Gray64:
+            return bytes_to_gray64Pixel;
         default:
             always_assert(NOT_IMPLEMENTED_YET);
             break;
@@ -68,47 +73,61 @@ class Matrix
         data.resize(rowCount * colCount);
     }
 
+    ~Matrix()
+    {
+        for (size_t i = 0; i < data.size(); i++)
+        {
+            if (data[i])
+            {
+                delete data[i];
+            }
+        }
+    }
+
     // Get number of rows.
     uint rows() const { return this->rowCount; }
 
     // Get number of cols.
     uint cols() const { return this->colCount; }
 
-    inline T &at(const uint &row, const uint &col)
+    inline void set(const uint &row, const uint &col, BasePixel *px)
     {
-        return this->data[((row * this->colCount) + col)];
+        data[((row * colCount) + col)] = px;
     }
 
-    // Get constant element reference at row and col.
-    inline const T &at(const uint &row, const uint &col) const
+    inline const BasePixel *get(const uint &row, const uint &col)
     {
-        return this->data[((row * this->colCount) + col)];
+        return data[((row * colCount) + col)];
     }
 
-    // Get constant element reference at index.
-    inline const T &atIndex(const size_t &index) const
+    inline void set(const ulong index, BasePixel *px)
     {
-        return this->data[index];
+        data[index] = px;
     }
 
-    bool operator==(const Matrix<T> &other)
+    inline const BasePixel *get(const ulong index)
     {
-        return equals(other);
+        return data[index];
     }
 
-    bool operator!=(const Matrix<T> &other)
+    bool operator==(const Matrix &other)
     {
-        return (!equals(other));
+        return equals(other);
     }
 
-    void parse_bytes_to_image(std::vector<byte> &bytes, const PixelType pt)
+    bool operator!=(const Matrix &other)
     {
-        auto parse_pixel = get_pixel_parse_function(pt);
-        size_t count = rowCount * colCount;
-        ulong bufferIndex = 0;
-        for (size_t i = 0; i < count; i++)
-        {
-            data[i] = parse_pixel(bytes, bufferIndex);
-        }
+        return (!equals(other));
     }
+
+    // void parse_bytes_to_image(std::vector<byte> &bytes, const PixelType pt)
+    // {
+    //     auto parse_pixel = get_pixel_parse_function(pt);
+    //     size_t count = rowCount * colCount;
+    //     ulong bufferIndex = 0;
+    //     for (size_t i = 0; i < count; i++)
+    //     {
+    //         data[i] = parse_pixel(bytes, bufferIndex);
+    //     }
+    // }
 };
diff --git a/czi-format/czi-parser/image/pixel_parse_functions.h b/czi-format/czi-parser/image/pixel_parse_functions.h
index 5161de7f7e883f1124c9302184dcf8cf9eb29991..f69097d8117709fb189cc69dcc048ecf8b64e18d 100644
--- a/czi-format/czi-parser/image/pixel_parse_functions.h
+++ b/czi-format/czi-parser/image/pixel_parse_functions.h
@@ -3,9 +3,14 @@
 #include "../custom_types.h"
 #include "../czi_parts/pixel_type.h"
 
-Bgr24Pixel bytes_to_bgr24pixel(std::vector<byte> &bytes, ulong &index)
+Bgr24Pixel bytes_to_bgr24Pixel(std::vector<byte> &bytes, ulong &index)
 {
     Bgr24Pixel pixel = Bgr24Pixel(bytes[index], bytes[index + 1], bytes[index + 2]);
     index += 3;
     return pixel;
+}
+
+Gray64Pixel bytes_to_gray64Pixel(std::vector<byte> &bytes, ulong &index)
+{
+    always_assert(NOT_IMPLEMENTED_YET);
 }
\ No newline at end of file
diff --git a/czi-format/czi-parser/image/pixel_structures.h b/czi-format/czi-parser/image/pixel_structures.h
index b9a2bd0ab7487fbc3b6b7e94a3cd9fc40a95bb21..7890b3618f6f03bec091d469b569d88f99fd8eb3 100644
--- a/czi-format/czi-parser/image/pixel_structures.h
+++ b/czi-format/czi-parser/image/pixel_structures.h
@@ -3,7 +3,20 @@
 
 struct BasePixel
 {
+    // True if pixel wasn't initialized with value or if is destructed.
     bool isEmpty = true;
+
+    virtual ~BasePixel()
+    {
+        isEmpty = true;
+    }
+
+    // Virtually check equalness of two pixels.
+    virtual bool equals(const BasePixel *other) const
+    {
+        always_assert(false && "This should be overriden!");
+        return false;
+    };
 };
 
 struct Gray64Pixel : BasePixel
@@ -20,6 +33,17 @@ struct Gray64Pixel : BasePixel
     bool operator==(const Gray64Pixel &b) const { return (value == b.value); }
     bool operator!=(const Gray64Pixel &b) const { return (value != b.value); }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Gray64Pixel *casted = dynamic_cast<const Gray64Pixel *>(other);
+        if (casted)
+        {
+            return (casted->value == value);
+        }
+        return false;
+    }
+
     Gray64Pixel &operator+=(const Gray64Pixel &b)
     {
         value += b.value;
@@ -53,6 +77,17 @@ struct Gray32FloatPixel : BasePixel
     bool operator==(const Gray32FloatPixel &b) const { return (value == b.value); }
     bool operator!=(const Gray32FloatPixel &b) const { return (value != b.value); }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Gray32FloatPixel *casted = dynamic_cast<const Gray32FloatPixel *>(other);
+        if (casted)
+        {
+            return (casted->value == value);
+        }
+        return false;
+    }
+
     Gray32FloatPixel &operator-=(const Gray32FloatPixel &b)
     {
         value -= b.value;
@@ -86,6 +121,17 @@ struct Gray32Pixel : BasePixel
     bool operator==(const Gray32Pixel &b) const { return (value == b.value); }
     bool operator!=(const Gray32Pixel &b) const { return (value != b.value); }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Gray32Pixel *casted = dynamic_cast<const Gray32Pixel *>(other);
+        if (casted)
+        {
+            return (casted->value == value);
+        }
+        return false;
+    }
+
     Gray32Pixel &operator+=(const Gray32Pixel &b)
     {
         value += b.value;
@@ -119,6 +165,17 @@ struct Gray16Pixel : BasePixel
     bool operator==(const Gray16Pixel &b) const { return (value == b.value); }
     bool operator!=(const Gray16Pixel &b) const { return (value != b.value); }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Gray16Pixel *casted = dynamic_cast<const Gray16Pixel *>(other);
+        if (casted)
+        {
+            return (casted->value == value);
+        }
+        return false;
+    }
+
     Gray16Pixel &operator-=(const Gray16Pixel &b)
     {
         value -= b.value;
@@ -140,6 +197,17 @@ struct Gray8Pixel : BasePixel
     bool operator==(const Gray8Pixel &b) const { return (value == b.value); }
     bool operator!=(const Gray8Pixel &b) const { return (value != b.value); }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Gray8Pixel *casted = dynamic_cast<const Gray8Pixel *>(other);
+        if (casted)
+        {
+            return (casted->value == value);
+        }
+        return false;
+    }
+
     Gray8Pixel &operator-=(const Gray8Pixel &b)
     {
         value -= b.value;
@@ -160,6 +228,17 @@ struct Bgr96FloatPixel : BasePixel
     Bgr96FloatPixel operator-(const Bgr96FloatPixel &other) const { return Bgr96FloatPixel(b - other.b, g - other.g, r - other.r); }
     Bgr96FloatPixel operator*(const float &v) const { return Bgr96FloatPixel(b * v, g * v, r * v); }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Bgr96FloatPixel *casted = dynamic_cast<const Bgr96FloatPixel *>(other);
+        if (casted)
+        {
+            return ((b == casted->b) && (g == casted->g) && (r == casted->r));
+        }
+        return false;
+    }
+
     float operator*(const Bgr96FloatPixel &other) const
     {
         float result = (b * other.b) + (g * other.g) + (r * other.r);
@@ -206,6 +285,17 @@ struct Bgr48Pixel : BasePixel
     Bgr48Pixel operator-(const Bgr48Pixel &other) const { return Bgr48Pixel(b - other.b, g - other.g, r - other.r); }
     Bgr96FloatPixel operator*(const float &v) const { return Bgr96FloatPixel((float)b * v, (float)g * v, (float)r * v); }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Bgr48Pixel *casted = dynamic_cast<const Bgr48Pixel *>(other);
+        if (casted)
+        {
+            return ((b == casted->b) && (g == casted->g) && (r == casted->r));
+        }
+        return false;
+    }
+
     int operator*(const Bgr48Pixel &other) const
     {
         int result = (b * other.b) + (g * other.g) + (r * other.r);
@@ -236,6 +326,17 @@ struct Bgr24Pixel : BasePixel
     Bgr48Pixel operator+(const Bgr24Pixel &other) const { return Bgr48Pixel(b + other.b, g + other.g, r + other.r); }
     Bgr24Pixel operator-(const Bgr24Pixel &other) const { return Bgr24Pixel(b - other.b, g - other.g, r - other.r); }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Bgr24Pixel *casted = dynamic_cast<const Bgr24Pixel *>(other);
+        if (casted)
+        {
+            return ((b == casted->b) && (g == casted->g) && (r == casted->r));
+        }
+        return false;
+    }
+
     Bgr96FloatPixel operator*(const float &v) const
     {
         return Bgr96FloatPixel((float)b * v, (float)g * v, (float)r * v);
@@ -270,13 +371,21 @@ struct Bgra32Pixel : BasePixel
     Bgra32Pixel() {}
     Bgra32Pixel(byte blue, byte green, byte red, byte alpha) : b(blue), g(green), r(red), a(alpha) { isEmpty = false; }
 
-    //Bgra32Pixel operator+(const Bgra32Pixel &b) const;
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Bgra32Pixel *casted = dynamic_cast<const Bgra32Pixel *>(other);
+        if (casted)
+        {
+            return ((b == casted->b) && (g == casted->g) && (r == casted->r) && (a == casted->a));
+        }
+        return false;
+    }
+
     Bgra32Pixel operator-(const Bgra32Pixel &other) const { return Bgra32Pixel(b - other.b, g - other.g, r - other.r, a - other.a); }
-    //Bgra32Pixel operator*(const float &other) const;
     int operator*(const Bgra32Pixel &other) const { return ((b * other.b) + (g * other.g) + (r * other.a) + (a * other.a)); }
     bool operator==(const Bgra32Pixel &other) const { return ((b == other.b) && (g == other.g) && (r == other.r) && (a == other.a)); }
     bool operator!=(const Bgra32Pixel &other) const { return !((b == other.b) && (g == other.g) && (r == other.r) && (a == other.a)); }
-    //Bgra32Pixel &operator+=(const Bgra32Pixel &other);
     Bgra32Pixel &operator-=(const Bgra32Pixel &other)
     {
         b -= other.b;
@@ -319,6 +428,17 @@ struct Gray64ComplexFloatPixel : BasePixel
     bool operator==(const Gray64ComplexFloatPixel &b) const { return ((imaginaryPart == b.imaginaryPart) && (realPart == b.realPart)); }
     bool operator!=(const Gray64ComplexFloatPixel &b) const { return !((imaginaryPart == b.imaginaryPart) && (realPart == b.realPart)); }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Gray64ComplexFloatPixel *casted = dynamic_cast<const Gray64ComplexFloatPixel *>(other);
+        if (casted)
+        {
+            return ((imaginaryPart == casted->imaginaryPart) && (realPart == casted->realPart));
+        }
+        return false;
+    }
+
     Gray64ComplexFloatPixel &operator+=(const Gray64ComplexFloatPixel &b)
     {
         imaginaryPart += b.imaginaryPart;
@@ -398,6 +518,22 @@ struct Bgr192ComplexFloatPixel : BasePixel
                 (rRealPart * b.rRealPart));
     }
 
+    bool equals(const BasePixel *other) const override
+    {
+        // We can check if other is of type Gray64Pixel by dynamic_cast, if it returns 0 it is not.
+        const Bgr192ComplexFloatPixel *casted = dynamic_cast<const Bgr192ComplexFloatPixel *>(other);
+        if (casted)
+        {
+            return ((casted->bImaginaryPart == casted->bImaginaryPart) &&
+                    (bRealPart == casted->bRealPart) &&
+                    (gImaginaryPart == casted->gImaginaryPart) &&
+                    (gRealPart == casted->gRealPart) &&
+                    (rImaginaryPart == casted->rImaginaryPart) &&
+                    (rRealPart == casted->rRealPart));
+        }
+        return false;
+    }
+
     bool operator==(const Bgr192ComplexFloatPixel &b) const
     {
         return ((bImaginaryPart == b.bImaginaryPart) &&
diff --git a/czi-format/czi-parser/main.cpp b/czi-format/czi-parser/main.cpp
index 6a15ac6dd6429e2d734f94e3567dbf639081d2b1..e2e5188ed9cddf340b31b8fb80c33f0d9b884b17 100644
--- a/czi-format/czi-parser/main.cpp
+++ b/czi-format/czi-parser/main.cpp
@@ -4,17 +4,6 @@
 
 int main(int argc, char **argv)
 {
-    Matrix<Gray8Pixel> bwImage = Matrix<Gray8Pixel>(5, 5);
-    bwImage.at(0, 0) = Gray8Pixel(100);
-    bwImage.at(0, 1) = Gray8Pixel(100);
-    bwImage.at(0, 2) = Gray8Pixel(100);
-    Matrix<Gray8Pixel> bwImage2 = Matrix<Gray8Pixel>(5, 5);
-    bwImage2.at(0, 0) = Gray8Pixel(100);
-    bwImage2.at(0, 1) = Gray8Pixel(101);
-    bwImage2.at(0, 2) = Gray8Pixel(100);
-    bool eq = bwImage == bwImage2;
-    printf("Matrices are %s\n", eq ? "equal" : "not equal");
-
     std::string cziFile;
     // #if DEBUG
     //     printf("***DEBUG MODE***\n");
@@ -58,5 +47,6 @@ int main(int argc, char **argv)
         parseResult.extract_images(dumpName);
 
     printf("Finished.\n");
+
     return 0;
 }
\ No newline at end of file