PNG: Fix check for too few axes, memory, orientations

core/file/png.cpp

@@ -183,6 +183,7 @@ namespace MR
           bit_depth (0),
           filename (filename),
           data_type (H.datatype()),
+          multiplier (1.0),
           outfile (NULL)
       {
         if (Path::exists (filename) && !App::overwrite_files)
@@ -196,6 +197,8 @@ namespace MR
           throw Exception ("Unable to set jump buffer for PNG structure for image \"" + filename + "\"");
         }
         outfile = fopen (filename.c_str(), "wb");
+        if (!outfile)
+          throw Exception ("Unable to open PNG file for writing for image \"" + filename + "\": " + strerror (errno));
         png_init_io (png_ptr, outfile);
         png_set_compression_level (png_ptr, Z_DEFAULT_COMPRESSION);
         switch (H.ndim()) {
@@ -231,17 +234,23 @@ namespace MR
             png_destroy_write_struct (&png_ptr, &info_ptr);
             throw Exception ("Undefined data type in image \"" + H.name() + "\" for PNG writer");
           case DataType::Bit:
-            bit_depth = 1;
+            assert (false);
             break;
           case DataType::UInt8:
+            bit_depth = 8;
+            break;
           case DataType::Float32:
             bit_depth = 8;
+            multiplier = std::numeric_limits<uint8_t>::max(); break;
             break;
           case DataType::UInt16:
           case DataType::UInt32:
           case DataType::UInt64:
+            bit_depth = 16;
+            break;
           case DataType::Float64:
             bit_depth = 16;
+            multiplier = std::numeric_limits<uint16_t>::max(); break;
             break;
         }
         // Detect cases where one axis has a size of 1, and hence represents the image plane
@@ -323,10 +332,11 @@ namespace MR
             row_pointers[row] = to_write + row * row_bytes;
           png_write_image (png_ptr, row_pointers);
           png_write_end (png_ptr, info_ptr);
+          delete[] row_pointers;
         };
 
 
-        if (bit_depth == 1 || data_type == DataType::UInt8 || data_type == DataType::UInt16BE) {
+        if (data_type == DataType::UInt8 || data_type == DataType::UInt16BE) {
           finish (data);
         } else {
           uint8_t scratch[row_bytes * height];

core/file/png.h

@@ -83,6 +83,7 @@ namespace MR
           int color_type, bit_depth;
           std::string filename;
           DataType data_type;
+          default_type multiplier;
           FILE* outfile;
 
           static void error_handler (png_struct_def*, const char*);
@@ -100,16 +101,8 @@ namespace MR
         std::function<default_type(const void*,size_t,default_type,default_type)> fetch_func;
         std::function<void(default_type,void*,size_t,default_type,default_type)> store_func;
         __set_fetch_store_functions<default_type> (fetch_func, store_func, data_type);
-        default_type multiplier = 1.0;
-        switch (data_type() & DataType::Type) {
-          case DataType::Float32: multiplier = std::numeric_limits<uint8_t>::max(); break;
-          case DataType::Float64: multiplier = std::numeric_limits<uint16_t>::max(); break;
-        }
-        for (size_t i = 0; i != num_elements; ++i) {
-          Raw::store_BE<T> (std::min (default_type(std::numeric_limits<T>::max()), std::max (0.0, std::round(multiplier * fetch_func (in_ptr, 0, 0.0, 1.0)))), out_ptr);
-          in_ptr += data_type.bytes();
-          out_ptr += sizeof(T);
-        }
+        for (size_t i = 0; i != num_elements; ++i)
+          Raw::store_BE<T> (std::min (default_type(std::numeric_limits<T>::max()), std::max (0.0, std::round(multiplier * fetch_func (in_ptr, i, 0.0, 1.0)))), out_ptr, i);
       };
 
 

core/formats/png.cpp

@@ -70,16 +70,16 @@ namespace MR
       }
 
       H.size(0) = png.get_width();
-      H.stride(0) = -3;
+      H.stride(0) = -2;
 
       H.size(1) = png.get_height();
-      H.stride(1) = -4;
+      H.stride(1) = -3;
 
       H.size(2) = 1;
-      H.stride(2) = 1;
+      H.stride(2) = 4;
 
       if (H.ndim() == 4)
-        H.stride(3) = 2;
+        H.stride(3) = 1;
 
       H.spacing (0) = H.spacing (1) = H.spacing (2) = 1.0;
       H.transform().setIdentity();
@@ -132,14 +132,14 @@ namespace MR
       if (H.ndim() == 4 && H.size(3) > 4)
         throw Exception ("A 4D image written to PNG must have between one and four volumes (requested: " + str(H.size(3)) + ")");
 
-      // TODO After looping over axes via square-bracket notation,
+      // After looping over axes via square-bracket notation,
       //   there needs to be at least two axes with size greater than one
-      size_t unity_axes = 0;
+      size_t nonunity_axes = 0;
       for (size_t axis = 0; axis != H.ndim(); ++axis) {
-        if (H.size (axis) == 1)
-          ++unity_axes;
+        if (H.size (axis) > 1)
+          ++nonunity_axes;
       }
-      if (unity_axes - (H.ndim() - num_axes) < 2)
+      if (nonunity_axes - (H.ndim() - num_axes) < 2)
         throw Exception ("Too few (non-unity) image axes to support PNG export");
 
       // For 4D images, can support:
@@ -149,7 +149,8 @@ namespace MR
       // - 4 volumes (save as RGBA)
       // This needs to be compatible with NameParser used in Header::create():
       //   "num_axes" subtracts from H.ndim() however many instances of [] there are
-      size_t width_axis = 0, axis_to_zero = 3;
+      // size_t width_axis = 0;
+      size_t axis_to_zero = 3;
       if (H.ndim() - num_axes > 1)
         throw Exception ("Cannot nominate more than one axis using square-bracket notation for PNG format");
       switch (num_axes) {
@@ -170,7 +171,7 @@ namespace MR
             axis_to_zero = 1;
           } else if (H.size(0) == 1) {
             axis_to_zero = 0;
-            width_axis = 1;
+            //width_axis = 1;
           } else {
             // If image is 3D, and all three axes have size greater than one, and we
             //   haven't used the square-bracket notation, we can't export genuine 3D data
@@ -192,8 +193,8 @@ namespace MR
           }
           if (axis < 0)
             throw Exception ("Cannot export 4D image to PNG format if all three spatial axes have size greater than 1 and square-bracket notation is not used");
-          if (!axis_to_zero)
-            width_axis = 1;
+          // if (!axis_to_zero)
+          //   width_axis = 1;
           break;
         default:
           throw Exception ("Cannot generate PNG file(s) from image with more than 4 axes");
@@ -209,7 +210,7 @@ namespace MR
       H.stride(1) = -3;
       H.spacing(1) = 1.0;
       if (H.ndim() > 2) {
-        H.stride(2) = 4;
+        H.stride(2) = -4;
         H.spacing(2) = 1.0;
       }
       if (H.ndim() > 3) {
@@ -223,9 +224,10 @@ namespace MR
 
       H.transform().setIdentity();
 
-      if (H.datatype() == DataType::Bit && H.size (width_axis) % 8) {
-        WARN ("Cannot write bitwise PNG image with width not a factor of 8; will instead write with 8-bit depth");
+      if (H.datatype() == DataType::Bit) {
+        WARN ("Cannot write bitwise PNG images; will instead write with 8-bit depth");
         H.datatype() = DataType::UInt8;
+        H.intensity_scale() = 1.0 / 255.0;
       }
 
       return true;

core/header.cpp

@@ -785,12 +785,13 @@ namespace MR
     Header result (headers[0]);
 
     if (axis_to_concat >= result.ndim()) {
+      Stride::symbolise (result);
       result.ndim() = axis_to_concat + 1;
       result.size(axis_to_concat) = 1;
+      result.stride(axis_to_concat) = axis_to_concat+1;
+      Stride::actualise (result);
     }
 
-    result.stride (axis_to_concat) = result.ndim()+1;
-
     for (size_t axis = 0; axis != result.ndim(); ++axis) {
       if (axis != axis_to_concat && result.size (axis) <= 1) {
         for (const auto& H : headers) {

core/image_io/png.cpp

@@ -29,17 +29,14 @@ namespace MR
 
     void PNG::load (const Header& header, size_t)
     {
-      segsize = header.datatype().bytes() * voxel_count (header) * files.size();
+      DEBUG (std::string("loading PNG image") + (files.size() > 1 ? "s" : "") + " \"" + header.name() + "\"");
       addresses.resize (1);
-      addresses[0].reset (new uint8_t [segsize]);
+      segsize = (header.datatype().bits() * voxel_count (header) + 7) / 8;
+      addresses[0].reset (new uint8_t[segsize]);
       if (is_new) {
         memset (addresses[0].get(), 0x00, segsize);
-        DEBUG ("allocated memory for PNG image \"" + header.name() + "\"");
       } else {
-        DEBUG (std::string("loading PNG image") + (files.size() > 1 ? "s" : "") + " \"" + header.name() + "\"");
-        size_t slice_bytes = (header.datatype().bits() * header.size(0) * header.size(1) + 7) / 8;
-        if (header.ndim() == 4)
-          slice_bytes *= header.size (3);
+        const size_t slice_bytes = (header.datatype().bits() * header.size(0) * header.size(1) * (header.ndim() == 4 ? header.size(3) : 1) + 7) / 8;
         for (size_t i = 0; i != files.size(); ++i) {
           File::PNG::Reader png (files[i].name);
           if (png.get_width() != header.size(0) ||
@@ -59,19 +56,15 @@ namespace MR
 
     void PNG::unload (const Header& header)
     {
-      if (addresses.size()) {
-        if (writable) {
-          size_t slice_bytes = (header.datatype().bits() * header.size(0) * header.size(1) + 7) / 8;
-          if (header.ndim() == 4)
-            slice_bytes *= header.size (3);
-          for (size_t i = 0; i != files.size(); i++) {
-            File::PNG::Writer png (header, files[i].name);
-            png.save (addresses[0].get() + (i * slice_bytes));
-          }
+      assert (addresses.size() == 1);
+      if (writable) {
+        const size_t slice_bytes = (header.datatype().bits() * header.size(0) * header.size(1) * (header.ndim() == 4 ? header.size(3) : 1) + 7) / 8;
+        for (size_t i = 0; i != files.size(); i++) {
+          File::PNG::Writer png (header, files[i].name);
+          png.save (addresses[0].get() + (i * slice_bytes));
         }
-        DEBUG ("deleting buffer for PNG image \"" + header.name() + "\"...");
-        addresses[0].release();
       }
+      delete[] addresses[0].release();
     }
 
   }

testing/binaries/tests/mrconvert

@@ -8,12 +8,22 @@ mrconvert mrconvert/in.mif tmp.nii  && testing_diff_image tmp.nii mrconvert/in.m
 mrconvert mrconvert/in.mif -datatype float32 tmp.nii.gz -force && testing_diff_image tmp.nii.gz mrconvert/in.mif
 mrconvert mrconvert/in.mif -strides 3,2,1 tmp.mgh  && testing_diff_image tmp.mgh mrconvert/in.mif
 mrconvert mrconvert/in.mif -strides 1,3,2 -datatype int16 tmp.mgz  && testing_diff_image tmp.mgz mrconvert/in.mif
-mrconvert mrconvert/in.mif tmp-[].png  && echo -e "1 0 0 0\n0 1 0 0\n0 0 1 0\n" > tmp.txt  && testing_diff_image tmp-[].png $(mrcalc mrconvert/in.mif 0 -max - | mrtransform - -replace tmp.txt -)
-mrconvert unit_warp.mif tmp-[].png -datatype uint8 -force  && echo -e "1 0 0 0\n0 1 0 0\n0 0 1 0\n" > tmp.txt  && testing_diff_image tmp-[].png $(mrcalc unit_warp.mif 0 -max -round - | mrtransform - -replace tmp.txt - | mrconvert - -vox 1,1,1 -)
+rm -f tmp-*.png && mrconvert mrconvert/in.mif tmp-[].png && testing_diff_image tmp-[].png $(mrcalc mrconvert/in.mif 0 -max - | mrtransform - -replace identity.txt -)
+rm -f tmp-*.png && mrconvert unit_warp.mif tmp-[].png -datatype uint8 && testing_diff_image tmp-[].png $(mrcalc unit_warp.mif 0 -max -round - | mrtransform - -replace identity.txt - | mrconvert - -vox 1,1,1 -)
 mrconvert dwi.mif tmp-[].mif -force && testing_diff_image dwi.mif tmp-[].mif
 mrconvert dwi.mif -coord 3 0:2:end tmp1.mif -force && mrconvert tmp-[0:2:66].mif tmp2.mif && testing_diff_header -keyval tmp1.mif tmp2.mif && testing_diff_image tmp1.mif tmp2.mif
 mrconvert mrcat/voxel[].mih - | testing_diff_header -keyval - mrcat/all_axis0.mif
 mrconvert mrcat/all_axis3.mif tmp-[].mif -force && testing_diff_header -keyval tmp-0.mif mrcat/voxel1.mih && testing_diff_header -keyval tmp-1.mif mrcat/voxel2.mih && testing_diff_header -keyval tmp-2.mif mrcat/voxel3.mih && testing_diff_header -keyval tmp-3.mif mrcat/voxel4.mih && testing_diff_header -keyval tmp-4.mif mrcat/voxel5.mih && testing_diff_header -keyval tmp-5.mif mrcat/voxel6.mih
 mrconvert dwi.mif tmp-[]-[].mif -force && testing_diff_image dwi.mif tmp-[]-[].mif
 mrconvert dwi.mif -coord 3 1:2:end -axes 0:2,-1,3 - | testing_diff_image - mrconvert/dwi_select_axes.mif
-
+mrinfo template.mif.gz -transform > tmp.txt && mrconvert template.mif.gz tmp[].png -force && mrconvert tmp[].png -vox 2.5 - | mrtransform - -replace tmp.txt - | mrcalc - 255 -div - | testing_diff_image - $(mrcalc template.mif.gz 1.0 -min -) -abs 0.002
+rm -f tmpaxial*.png && mrconvert template.mif.gz -coord 2 9,19,29,39,49 tmpaxial[].png && testing_diff_image tmpaxial[].png mrconvert/pngaxial[].png
+mrconvert mrconvert/pngaxial[].png - | testing_diff_image - mrconvert/pngaxial.mif.gz
+rm -f tmpcoronal*.png && mrconvert template.mif.gz -coord 1 17,32,47,62,77 -axes 0,2,1 tmpcoronal[].png && testing_diff_image tmpcoronal[].png mrconvert/pngcoronal[].png
+mrconvert mrconvert/pngcoronal[].png - | testing_diff_image - $(mrconvert mrconvert/pngcoronal.mif.gz -axes 0,2,1 -)
+rm -f tmpsagittal*.png && mrconvert template.mif.gz -coord 0 27,47,67 -axes 1,2,0 tmpsagittal[].png && testing_diff_image tmpsagittal[].png mrconvert/pngsagittal[].png
+mrconvert mrconvert/pngsagittal[].png - | testing_diff_image - $(mrconvert mrconvert/pngsagittal.mif.gz -axes 1,2,0 -)
+rm -f tmpmask*.png && mrconvert mask.mif tmpmask[].png && testing_diff_image tmpmask[].png mrconvert/pngmask[].png
+mrconvert mrconvert/pngmask[].png - | testing_diff_image - mrconvert/pngmask.mif.gz
+rm -f tmptissues*.png && mrconvert dwi2fod/msmt/tissues.mif tmprgb[].png && testing_diff_image tmprgb[].png mrconvert/pngrgb[].png
+mrconvert mrconvert/pngrgb[].png - | testing_diff_image - $(mrconvert dwi2fod/msmt/tissues.mif -vox 1,1,1 - | mrtransform - -replace identity.txt - | mrcalc - 255 -mult -round 255 -min -)