Save 2D double arrays as 3-channel png files (RGB)

I have a 2D data matrix (doubles) that I need to split in smaller frames and save as both .txt and .png files. The following is the routine that I wrote, and it has been working well until I applied it on a particular dataset (.txt input file is too big to be uploaded here).

After I load the input text file, I split them in operlapping frames, and save the output frames as .txt and .png. To save it as RGB image, I am using colormap extracted via 'imagesc', and a matlab function called double2rgb (downloaded from Matlab file exchange, link below). For this particular input file, some of the saved png images do not match the original imagesc patterns (please see the PNG files for the full frame (imagesc) as well as a screenshot of some of the saved split frames (1-36) that shows the issue).

If you look at the full frame, the upper left (light green bcakground, which are not constant values but very small values close to zero) - my guess is that this is the issue, because that is where the saved split frame images do not match the corresponding loacation pattern in the fullframe display (please look at the images attached here for comparison) - notice the images 'testout_1.png' through 'testout_19.png' in the screenshot. I appreciate any help/comment on this regard. Thanks!

https://www.mathworks.com/matlabcentral/fileexchange/30264-double2rgb

data = 'test.txt';    % Set input data
outdir_250 = '092423-split_frames_250_test';    % Set output directory
% Split frame (save as PNG files)
splitFrame(data,250,250,50,'testout',"png",outdir_250)
%% ====================================================================

function splitFrame(img,img_width,img_height,overlap,out,fext,outdir)

% Function input parameters:
% img = input file/image
% img_width = split frame/image width
% img_height = split frame/image height
% overlap = number of overlap data-points/pixel for sliding window 
% out = output filename prefix
% fext = output filename extention
% outdir = directory/folder where the split frames/images are saved 
%          (this will create a directory if it does not exist)
% 
% example:
% splitFrame(data,250,250,50,'South_Atlantic',"txt",outdir_south_atlantic_250)

img = load (img);  % Load data

x = size(img,1);  % Image x-dimension
y = size(img,2);  % Image y-dimension

imagesc(img)
cmap=colormap;    % get current colormap

% Create output directory, if it does not exist
if not(isfolder(outdir))
    mkdir(outdir)
end

if fext=="png"
    im = zeros(img_width,img_height,3);  % Initialize an emty 3D-array for split images (PNG)
else
    im = zeros(img_width,img_height);  % Initialize an emty 2D-array for split images (text)  
end

k = 1;  % Set a counter for naming the split images

for i = 1:overlap:x
    if (i+img_width) > x
        break
    end
    for j = 1:overlap:y

        if (j+img_height) > y
            break
        end

        im = img(i:i+img_width-1,j:j+img_height-1,:);
        filename = out + "_" + k + "." + fext;

        if fext=="png"
            imwrite(double2rgb(im,cmap),fullfile(outdir,filename));  % Save rgb split image/matrix          
        else
            save(fullfile(outdir,filename), 'im', '-ascii')
        end
        
        clear im
        k = k+1;

    end    
end

end

Answer 1

[DGM]

The file double2rgb() does not do the same quantization as pcolor()/imagesc() will with CDataScaling = 'scaled'.

Consider this example comparing double2rgb() against imagesc() and MIMT gray2pcolor().

% a color map
map = parula(16);

% an image on an arbitrary scale
A = lingrad([80 500],[0 0; 0 1],[0; 1],'double'); % unit-scale [0 1]
A = A*12 - 3.5; % [-3.5 8.5]

% quantize and apply the colormap
B = gray2pcolor(A,map,'default'); % same quantization as imquantize()
C = gray2pcolor(A,map,'cdscale'); % same quantization as imagesc()/pcolor()
D = double2rgb(double(A),map); % this will fail if A isn't already float

% capture the output of imagesc() just for comparison
imagesc(A)
colormap(map)
set(gca,'TickLength',[0 0]);
screenshot = frame2im(getframe(gca));
screenshot = imresize(screenshot,size(A),'nearest');
screenshot = im2double(screenshot);
clf

% compare the three tools against the screenshot
outpict = [B; C; D; screenshot];
imshow2(outpict)

Regarding double2rgb(), I would go so far as to say that this is less of a "difference" and more of a simple OBOE. I dunno; maybe objectives were different than what I expect.

The 'cdscale' mode of gray2pcolor()/uniquant() is intended to emulate imagesc(). While it's close, the operations do differ, so there is always room for slight differences in the way things get rounded.

gray2pcolor() is part of MIMT, though I have posted compact versions on the forum:

https://www.mathworks.com/matlabcentral/answers/1946888-how-to-remove-white-borders-space-from-an-rgb-image#answer_1216073

The second part of the problem is that you're relying on the input limits to be selected implicitly. The quantization process involves normalizing the data, typically with respect to the data extrema. If you process the image in chunks, the extrema of each chunk will vary, so the mapping will vary too. If you want the chunks to appear like the original whole image, you need to process them with respect to the extrema of the whole image. Both gray2pcolor() and double2rgb() allow these to be specified explicitly.

% a color map
map = parula(16);

% an image on an arbitrary scale
A = lingrad([80 500],[0 0; 0 1],[0; 1],'double'); % unit-scale [0 1]
A = A*12 - 3.5; % [-3.5 8.5]

% cut out only the left-hand half of the image
Alh = A(:,1:250);

% quantize and apply the colormap
limits = imrange(A); % get the extrema of the entire image
Clocal = gray2pcolor(Alh,map,'cdscale'); % use only the extrema of the LH half (keeps changing!)
Cglobal = gray2pcolor(Alh,map,limits,'cdscale'); % use the extrema of the entire image

% capture the output of imagesc() just for comparison
imagesc(Alh)
colormap(map)
set(gca,'TickLength',[0 0]);
caxis(limits); % specify this explicitly just the same
screenshot = frame2im(getframe(gca));
screenshot = imresize(screenshot,size(Alh),'nearest');
screenshot = im2double(screenshot);
clf

% compare the results
outpict = [Clocal; Cglobal; screenshot];
imshow2(outpict)

Alternatively, you might consider converting the entire image to RGB first and then detiling it. I don't know if that works for your needs.