v4.0.0: Rectangular arrays for image data (~31% faster), fix read-only NullRef

Image data now uses rectangular multi-dimensional arrays (float[,], float[,,])
instead of jagged arrays, eliminating ~12,500 heap allocations and enabling
zero-copy SIMD reads on .NET 10+ via MemoryMarshal.GetArrayDataReference.
299MB test image: 285ms -> 197ms.

Also fixes NullReferenceException when closing/flushing read-only files
(BufferedDataStream._out null guard) and incorrect FileMode.OpenOrCreate
for read-only BufferedFile access.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: sebgod

.gitignore

@@ -1,6 +1,11 @@
 ## Ignore Visual Studio temporary files, build results, and
 ## files generated by popular Visual Studio add-ons.
 
+# Large FITS test files (download separately)
+tests/CSharpFITS.Test/testdocs/LDN1089_singleFrame.fits
+*.local.json
+*.lock.json
+
 *.cobertura.xml
 
 # User-specific files

CSharpFITS/CSharpFITS.csproj

@@ -11,13 +11,13 @@
     <LangVersion>9.0</LangVersion>
     <PackageReadmeFile>README.md</PackageReadmeFile>
     <PackageLicenseFile>license.txt</PackageLicenseFile>
-    <AssemblyVersion>3.2.0.0</AssemblyVersion>
-    <FileVersion>3.2.0.0</FileVersion>
+    <AssemblyVersion>4.0.0.0</AssemblyVersion>
+    <FileVersion>4.0.0.0</FileVersion>
     <PackageTags>Astronomy,FITS,Image</PackageTags>
     <RepositoryUrl>https://github.com/SharpAstro/FITS.Lib</RepositoryUrl>
     <PackageProjectUrl>https://github.com/SharpAstro/FITS.Lib</PackageProjectUrl>
     <Description>Supports reading and writing FITS files</Description>
-    <Version>3.2.0</Version>
+    <Version>4.0.0</Version>
     <Authors>Virtual Observatory - India</Authors>
     <PackageId>FITS.Lib</PackageId>
     <RootNamespace>nom.tam</RootNamespace>

CSharpFITS/fits/ImageData.cs

@@ -338,7 +338,7 @@ public override void Read(ArrayDataIO i)
 			}
 			else
 			{
-				dataArray = ArrayFuncs.NewInstance(dataDescription.type, dataDescription.dims);
+				dataArray = ArrayFuncs.NewRectangularInstance(dataDescription.type, dataDescription.dims);
 				try
 				{
 					i.ReadArray(dataArray);
@@ -395,7 +395,7 @@ public override void Write(ArrayDataIO o)
 				else if (dataArray == null && dataDescription != null)
 				{
 					// Need to create an array to match a specified header.
-					dataArray = ArrayFuncs.NewInstance(dataDescription.type, dataDescription.dims);
+					dataArray = ArrayFuncs.NewRectangularInstance(dataDescription.type, dataDescription.dims);
 				}
 				else
 				{

CSharpFITS/image/ImageTiler.cs

@@ -37,7 +37,7 @@ public virtual Object CompleteImage
                     throw new IOException("Attempt to read from null file");
                 }
                 long currentOffset = f.Position;
-                Object o = ArrayFuncs.NewInstance(base_Renamed, dims);
+                Object o = ArrayFuncs.NewRectangularInstance(base_Renamed, dims);
                 f.Seek(fileOffset, SeekOrigin.Begin);
                 f.ReadArray(o);
                 f.Seek(currentOffset, SeekOrigin.Begin);
@@ -232,19 +232,16 @@ protected internal virtual void FillTile(Array data, Array o, int[] dims, int[]
         protected internal virtual void FillMemData(Array data, int[] posits, int length,
                                                 Array output, int outputOffset, int dim)
         {
-            // FIX THIS n-D crap
-            //if(data is Object[])
             if (ArrayFuncs.CountDimensions(data) > 1)
             {
                 if (ArrayFuncs.IsArrayOfArrays(data))
                 {
-                    //Object[] xo = (Object[]) data;
-                    //FillMemData((Array)xo[posits[dim]], posits, length, output, outputOffset, dim + 1);
                     FillMemData((Array)data.GetValue(posits[dim]), posits, length, output, outputOffset, dim + 1);
                 }
                 else
                 {
-                    throw new Exception("Called FillMemData with multi-dimensional array.");
+                    // Rectangular multi-dimensional array: compute linear offset and use Buffer.BlockCopy
+                    FillMemDataRectangular(data, posits, length, output, outputOffset);
                 }
             }
             else
@@ -269,6 +266,43 @@ protected internal virtual void FillMemData(Array data, int[] posits, int length
             }
         }
 
+        /// <summary>Fill tile data from a rectangular multi-dimensional array.</summary>
+        private void FillMemDataRectangular(Array data, int[] posits, int length,
+                                            Array output, int outputOffset)
+        {
+            int lastDim = dims.Length - 1;
+            int startFrom = posits[lastDim];
+            int startTo = outputOffset;
+            int copyLength = length;
+
+            if (posits[lastDim] < 0)
+            {
+                startFrom -= posits[lastDim];
+                startTo -= posits[lastDim];
+                copyLength += posits[lastDim];
+            }
+            if (posits[lastDim] + length > dims[lastDim])
+            {
+                copyLength -= (posits[lastDim] + length - dims[lastDim]);
+            }
+
+            if (copyLength <= 0) return;
+
+            // Compute linear offset into the rectangular array
+            int linearOffset = 0;
+            int stride = 1;
+            for (int d = lastDim; d >= 0; d--)
+            {
+                linearOffset += posits[d] < 0 ? 0 : posits[d] * stride;
+                stride *= dims[d];
+            }
+            // Adjust for the startFrom correction
+            linearOffset += (startFrom - (posits[lastDim] < 0 ? 0 : posits[lastDim]));
+
+            int elementSize = System.Runtime.InteropServices.Marshal.SizeOf(base_Renamed);
+            Buffer.BlockCopy(data, linearOffset * elementSize, output, startTo * elementSize, copyLength * elementSize);
+        }
+
         /// <summary>File a tile segment from a file.</summary>
         /// <param name="output">The output tile.</param>
         /// <param name="delta">The offset from the beginning of the image in bytes.</param>

CSharpFITS/util/ArrayFuncs.cs

@@ -142,27 +142,31 @@ public static int ComputeSize(Object o)
                 if (t.IsArray)
                 {
                     int size = 0;
+                    var array = (Array)o;
 
-                    if (t.HasElementType && t.GetElementType().IsPrimitive)
+                    // Rectangular multi-dimensional array: total elements * element size
+                    if (array.Rank > 1 && t.HasElementType && t.GetElementType().IsPrimitive)
+                    {
+                        size = array.Length * (int)sizes[t.GetElementType()];
+                    }
+                    else if (t.HasElementType && t.GetElementType().IsPrimitive)
                     {
-                        var array = (Array)o;
                         if (array.Length > 0)
                         {
-                            IEnumerator i = ((Array)o).GetEnumerator();
+                            IEnumerator i = array.GetEnumerator();
                             i.MoveNext();
 
                             size += array.Length * ComputeSize(i.Current);
                         }
                     }
                     else
                     {
-                        for (IEnumerator i = ((Array)o).GetEnumerator(); i.MoveNext();)
+                        for (IEnumerator i = array.GetEnumerator(); i.MoveNext();)
                         {
                             size += ComputeSize(i.Current);
                         }
                     }
 
-                    //return size;
                     result = size;
                 }
                 else if (t == typeof(String))
@@ -444,27 +448,23 @@ public static Type GetBaseClass(Object o)
 
             if (o.GetType().IsArray)
             {
-                // if 'o' is a jagged array of 2 or more than 2 dimensions
-                // with all dimensions of type premitive, 
-                // check for the type of last dimension.
-                if (o.GetType().ToString().StartsWith("System.Array"))
+                // Rectangular multi-dimensional array (e.g., float[,,])
+                if (((Array)o).Rank > 1)
                 {
-                    return GetBaseClass(((Array)o).GetValue(0));
+                    return o.GetType().GetElementType();
                 }
-                // if 'o' is a multidimension array of any dimension 
-                // or a jagged array of 2 or more than 2 dimension having
-                // first dimension of non-primitive type,
-                // return the type of 'o' directly without dimension brackets.
-                else
+                // Jagged array of arrays (e.g., Array[] of Array[] of float[])
+                if (o.GetType().ToString().StartsWith("System.Array"))
                 {
-                    return Type.GetType(o.ToString().Substring(0, o.ToString().IndexOf("[")));
+                    return GetBaseClass(((Array)o).GetValue(0));
                 }
+                // 1D or jagged typed array - strip all bracket notation to get base type
+                return Type.GetType(o.ToString().Substring(0, o.ToString().IndexOf("[")));
             }
             else
             {
                 return o.GetType();
             }
-
         }
 
         /// <summary>This routine returns the base class of an object with the dimension indication
@@ -582,19 +582,27 @@ public static Object Flatten(Object input)
         protected internal static int DoFlatten(Object input, Array output, int offset)
         {
             int result = 0;
+            var array = (Array)input;
+
+            // Rectangular multi-dimensional array: use Buffer.BlockCopy for efficiency
+            if (array.Rank > 1 && array.GetType().GetElementType().IsPrimitive)
+            {
+                int elementSize = (int)sizes[array.GetType().GetElementType()];
+                Buffer.BlockCopy(array, 0, output, offset * elementSize, array.Length * elementSize);
+                return array.Length;
+            }
 
             if (IsArrayOfArrays(input))
             {
-                int i = offset;
-                for (IEnumerator e = ((Array)input).GetEnumerator(); e.MoveNext();)
+                for (IEnumerator e = array.GetEnumerator(); e.MoveNext();)
                 {
                     Object a = e.Current;
                     result += DoFlatten(a, output, offset + result);
                 }
             }
             else if (input.GetType().IsArray)
             {
-                IEnumerator e = ((Array)input).GetEnumerator();
+                IEnumerator e = array.GetEnumerator();
                 for (int i = offset; e.MoveNext(); ++i)
                 {
                     output.SetValue(e.Current, i);
@@ -643,19 +651,13 @@ public static Array Curl(Array input, int[] dimens)
 
             for (IEnumerator i = input.GetEnumerator(); i.MoveNext() && NextIndex(index, dimens);)
             {
-                //NewInstance call creates a jagged array. So we cannot set the value using Array.SetValue
-                //as it works for multi-dimensional arrays and not for jagged
-                //newArray.SetValue(i.Current, index);
-
                 Array tarr = newArray;
                 for (int i2 = 0; i2 < index.Length - 1; i2++)
                 {
                     tarr = (Array)tarr.GetValue(index[i2]);
                 }
                 tarr.SetValue(i.Current, index[index.Length - 1]);
             }
-            //int offset = 0;
-            //DoCurl(input, newArray, dimens, offset);
             return newArray;
         }
         /// <summary>
@@ -791,22 +793,26 @@ public static Array NewInstance(Type cl, int[] dims, int dim)
 
                 return a;
             }
-            /*
-            Array o = Array.CreateInstance(cl, dims);
-            if(o == null)
+        }
+
+        /// <summary>Create a rectangular (multi-dimensional) array with a single allocation.
+        /// For 1D arrays, returns a simple typed array. For 2D+, returns e.g. float[,] or float[,,].
+        /// This is more efficient than jagged arrays for image data.</summary>
+        /// <param name="cl">The element type.</param>
+        /// <param name="dims">The dimensions.</param>
+        public static Array NewRectangularInstance(Type cl, int[] dims)
+        {
+            if (dims.Length == 0)
             {
-              System.String desc = cl + "[";
-              System.String comma = "";
-              for(int i = 0; i < dims.Length; i += 1)
-              {
-                desc += comma + dims[i];
-                comma = ",";
-              }
-              desc += "]";
-              throw new System.OutOfMemoryException("Unable to allocate array: " + desc);
+                dims = new int[] { 1 };
             }
-            return o;
-            */
+
+            if (dims.Length == 1)
+            {
+                return NewInstance(cl, dims[0]);
+            }
+
+            return Array.CreateInstance(cl, dims);
         }
 
 
@@ -869,6 +875,27 @@ public static bool ArrayEquals(Object x, Object y, double tolf, double told)
 
             if (xClass != yClass)
             {
+                // Allow comparing jagged vs rectangular arrays with the same base type and dimensions
+                if (x is Array xa && y is Array ya)
+                {
+                    Type xBase = GetBaseClass(x);
+                    Type yBase = GetBaseClass(y);
+                    int[] xDims = GetDimensions(x);
+                    int[] yDims = GetDimensions(y);
+                    if (xBase == yBase && xDims.Length == yDims.Length)
+                    {
+                        bool sameDims = true;
+                        for (int i = 0; i < xDims.Length; i++)
+                            if (xDims[i] != yDims[i]) { sameDims = false; break; }
+                        if (sameDims && xBase.IsPrimitive)
+                        {
+                            // Flatten both and compare bytes
+                            Object fx = Flatten(x);
+                            Object fy = Flatten(y);
+                            return ArrayEquals(fx, fy, tolf, told);
+                        }
+                    }
+                }
                 return false;
             }