The Unremarkable Adventures of an Electron: July 2015

Wednesday, July 29, 2015

Read and write animated gifs with MATLAB

UPDATE: These scripts have been vastly improved. Find the current versions here.

In the process of developing my own image mangling toolbox for Matlab, I had routine need for generating animated gif files from image sequences. It then follows that said .gif files might need to be read back again. Matlab's inbuilt functions imread() and imwrite() don't make this trivial. There are some suggestions floating around the internet, and there are some unanswered reports of unexpected behavior. Not everything worked simply, but I came up with my own ways.

First, writing the image was fairly simple. I had originally been using imagemagick to do the heavy lifting, but decided to go for a more direct route. As mentioned, the imagemagick method does seem to have better output, but it's quite a bit slower in my experience.

function gifwrite(inarray,filepath,delay,method)
%   GIFWRITE(INARRAY, FILEPATH, {DELAY}, {METHOD})
%       write image stack to an animated gif
%       
%   INARRAY: 4-D image array (rgb, uint8)
%   FILEPATH: full name and path of output animation
%   DELAY: frame delay in seconds (default = 0.05)
%   METHOD: animation method, 'native' or 'imagemagick' (default = 'native')
%       'imagemagick' may have better quality, but is much slower
 
if nargin<4;
    method='native';
end
 
if nargin<3;
    delay=0.05;
end
 
numframes=size(inarray,4);
 
if strcmpi(method,'native');
    disp('creating animation')
    for n=1:1:numframes;
        [imind,cm]=rgb2ind(inarray(:,:,:,n),256);
        if n==1;
            imwrite(imind,cm,filepath,'gif','DelayTime',delay,'Loopcount',inf);
        else
            imwrite(imind,cm,filepath,'gif','DelayTime',delay,'WriteMode','append');
        end
    end
else
    disp('creating frames')    
    for n=1:1:numframes;
        imwrite(inarray(:,:,:,n),sprintf('/dev/shm/%03dgifwritetemp.png',n),'png');
    end
     
    disp('creating animation')
    system(sprintf('convert -delay %d -loop 0 /dev/shm/*gifwritetemp.png %s',delay*100,filepath));
     
    disp('cleaning up')    
    system('rm /dev/shm/*gifwritetemp.png');
end
return

Reading the image back wasn't straightforward at all. The naive approach suggested by the documentation and posts online does not result in correct output. Imread(...'frames','all') returns only a single colormap corresponding to the global color table in the file. Any file containing multiple images with local color tables will turn into a pile of garbage.

function outpict=gifreadcrap(filepath)
%   GIFREADCRAP(FILEPATH)
%       reads all frames of an animated gif into a 4-D RGB image array
%       seems imread() cannot correctly read animated gifs
 
[images map]=imread(filepath, 'gif','Frames','all');
 
s=size(images);
numframes=s(4);
 
outpict=zeros([s(1:2) 3 numframes],'uint8');
for n=1:1:numframes;
    outpict(:,:,:,n)=ind2rgb8(images(:,:,:,n),map);
end
 
return


Original image written with gifwrite() and results as returned by gifreadcrap()

Using imread() to read single frames in the hope of perhaps getting correct color information produced the exact same results. Imfinfo() can be used to fetch file information, and certainly, it returns color tables which should correspond to the LCTs in the file. This is where things get ugly.

I spent some time with a hex editor and the file standards documentation to verify what I suspected was a bug. While imfinfo() returned the LCT data, they were all shifted by exactly one byte. The immediately adjacent bytes were being read correctly, though. Sounds like an OBOE to me!

At this point, I discover that yes, indeed it is a bug in versions R14-2012a. A patch exists, but for the sake of anyone who doesn't care, I decided to integrate both the native solution and my own imagemagick workaround.

Keep in mind that it may still be necessary to coalesce the animation before opening it, depending on what you want to do in Matlab. If you're trying to import an optimized gif, you can use the included option to coalesce the image. This optional mode is similar to the other optional modes in these two functions in that it requires external tools and assumes a linux environment. All temporary file operations utilize /dev/shm for marginal speed improvement. Altering this temporary path for your own environment should be trivial. Both functions should work in default modes on other systems, but I have no intention of testing that.

function outpict=gifread(filepath,method,coalesce)
%   GIFREAD(FILEPATH, {METHOD}, {COALESCE})
%       reads all frames of an animated gif into a 4-D RGB image array
%       
%   FILEPATH: full path and filename
%   METHOD: file read method, 'native' or 'imagemagick' (default = 'native')
%       'imagemagick' method is a workaround for bug 813126 present in
%       R14SP3-2012a versions.  Bug consists of an OBOE in reading LCT data.
%       A patch does exist for these versions:
%       https://www.mathworks.com/support/bugreports/813126
%   COALESCE: 0 or 1, Specifies whether to coalesce the image sequence prior to
%       importing.  Used when loading optimized gifs. Requires imagemagick.
%       (optional, default 0)
  
if nargin<3;
    coalesce=0;
end
if nargin<2;
    method='native';
end
  
  
if coalesce==1
    system(sprintf('convert %s -layers coalesce /dev/shm/gifreadcoalescetemp.gif',filepath));
    filepath='/dev/shm/gifreadcoalescetemp.gif';
end
  
if strcmpi(method,'native')
    % use imread() directly (requires patched imgifinfo.m)
    [images map]=imread(filepath, 'gif','Frames','all');
    infostruct=imfinfo(filepath);
  
    s=size(images);
    numframes=s(4);
  
    outpict=zeros([s(1:2) 3 numframes],'uint8');
    for n=1:1:numframes;
        LCT=infostruct(1,n).ColorTable;
        outpict(:,:,:,n)=ind2rgb8(images(:,:,:,n),LCT);
    end
else
    % split the gif using imagemagick instead
    system(sprintf('convert %s /dev/shm/%%03d_gifreadtemp.gif',filepath));
    [~,numframes]=system('ls -1 /dev/shm/*gifreadtemp.gif | wc -l');
  
    numframes=str2num(numframes);
    [image map]=imread('/dev/shm/000_gifreadtemp.gif', 'gif');
    s=size(image);
  
    outpict=zeros([s(1:2) 3 numframes],'uint8');
    for n=1:1:numframes;
        [image map]=imread(sprintf('/dev/shm/%03d_gifreadtemp.gif',n-1), 'gif');
        outpict(:,:,:,n)=ind2rgb8(image,map);
    end
  
    system('rm /dev/shm/*gifreadtemp.gif');
end
  
if coalesce==1
    system(sprintf('rm %s',filepath));
end
  
return

While all of this works, the native reading method does require either a patched copy of imgifinfo.m or Matlab version 2012b or later. Of course, once the existence of the bug was known, fixing things was simple. All the hours over the last day and a half that I spent digging for answers online and with a hex editor were to only rediscover something that was already known but hidden behind MathWorks' member login. It kind of pisses me off enough that finding the solution does not resolve my focus. What possible purpose does restricting access to bug reports serve?

If there's one thing I should have learned from my experiences with asking things of forums, it's to never ask forums. The internet is littered with the evidence of my inability to learn this simple lesson. This very blog was an angry reaction to the previous spectacularly infuriating experience. If I'm bound to ask questions of a silent screen -- if I'm bound to carve my own conclusions and place them on them complete on someone else's mantle in the meager hope that I can help the next person avoid my fate, then I might as well do it in a squalor of my own crafting, without the unrealistic expectations of interaction tugging at my attention.

Monday, July 20, 2015

Toward a useful vision aid

Ever since cataract surgery, I've been struggling with being able to see well. Part of the issue is a gradual increase in astigmatism since surgery (I just need to get my prescription updated), but the biggest hurdle is the loss of visual accommodation.

With the extraction of the eye's natural lens, goes the ability to adjust focus to accommodate for distance. While corrective lenses can bring things into focus for any fixed distance (e.g. reading glasses), a fixed correction can't work over a broad range of distances. As depth of field increases as a function of distance to subject, it stands to reason that most accommodation issues occur when working with nearby objects. An older person feeling the effects of presbyopia can probably say a thing or two about the gradual frustration of needing bifocals.

This isn't about age-related problems or old people things. This is about my dumb broke ass trying to find a means to get by with my pointless daily existence. Anyone else would probably just buy some magnifier lamps or something, but those cost double-digits kind of money. If you've been paying attention, you'll know that I'll end up cobbling shit out of $0.99 Ebay garbage and then sitting in the dark at midnight writing rambling stories to myself about an experience that even I can barely care about. That's what I did. That's what I'm doing. Let us continue.

Presbyopia is common though. Where do most people start? Bifocals and reading glasses simply offer an increased optical power (increased lens convexity) so that focusing on near objects is possible. While bifocals are convenient in that they're always at hand, their use requires a strong downward gaze which I can't effectively accomplish -- a consequence of an entirely separate ailment.

A drift of reading glasses and very uncomfortable Optivisor clone

Drug store reading glasses aren't meant to be used with prescription glasses, though it's ... possible. Prescription reading glasses are a simple addition to the SPH portion of the prescription (or specification of a +ADD). Of course, if you know your prescription and can use cheap readers to get an idea how much correction you need for a particular distance, you can do the math yourself and order any sort of odd double-reader bifocals or "computer glasses" you want.

For instance, let's say I wear plain single-vision glasses and I pick up a pair of +2.5 readers and slap them over the top. "Hooray!" I say as I can again read my own handwriting at 12". I can either add 2.5 to the SPH portion of both OD and OS lines on my prescription, or if it's appropriate, I can just use the +ADD entry on the order form. In this way, you can order your own readers or bifocals online without dealing with extra expense.

Say I have an existing bifocal prescription and I want to make computer glasses that can focus at 30", but without changing the power of the secondary lens I normally use for reading at 15". If I can determine that a +1.50 correction allows me to see at 30", I can just add that to the SPH section of the prescription and then subtract it from the existing +ADD which specified the original bifocals. Zenni Optical actually covers these sorts of prescription adjustments in their FAQ here and here.

Granted, buying glasses from Zenni beats paying $200 for glasses, but I'm not made of money. Not only that, but like I mentioned, depth of field is a function of distance. With glasses configured to focus at 20', you can focus at 200', but with glasses configured to focus at 2', you probably won't be able to see much at 20'. I went though the motions I describe above when I got my glasses for the computer. In fact, I collected the optimal focusing distance for a range of positive correction powers, as well as the minimum and maximum distance I could reasonably make out detailed edges in some text samples. It seems subjective at a glance, but when plotted, the trend starts to reveal a nice rational relationship.

With some use-cases such as the computer where distance is fixed and repeatable, it's simple to come up with a correction that's tailored to the task. Otherwise, shop work involving operations up close (checking pitch of a bolt with a thread gage) and at a distance (finding the bolt after you drop it) is hard to correct with a single prescription. At some point one has to accept that there is no single solution, and that carrying numerous solutions around is impractical.

So what's left? If I can't accommodate and I can't effectively use bifocals, I'm pretty much doomed to do half of my work in a blurry world unless I'm constantly swapping glasses or fumbling a magnifier. Some work patterns are harder to deal with than others. Automotive work requires focus between about 12" and 60", and the area of focus is typically changing regularly. On the other hand, lathe operations cover a smaller range of distances and most importantly, the critical areas of focus don't change as much. It's simpler to try coming up with task-oriented magnification tools for cases like this. ... simpler to try.

So I embark on a journey to come up with some positionable magnetic work magnifiers that I can stick down on the lathe or drill press or vise when I need them. There are lots of options out there already if you want to spend money. Certainly, I'd like a lighted magnifier on a Noga mag base, but that's not going to happen.

I found a MagniStitch magnifier that my late grandparents had purchased in the 80's. The lens is well-shaped and wide. The ball-jointed arm is much less useful than one might think. It tends to pop apart when positioning. I glued it to an old speaker magnet (because double-sided foam tape is bullshit). These can still be purchased for about $13-$20. I've considered buying others if perhaps I can CAD up a more appropriate articulated arm for it.

MagniStitch stuck on the drill press vise

The MagniStitch didn't work well for me on the lathe simply because the arm is too short to reach around the tool post or to reach around the cross slide to see the dial. I considered making my own from an inexpensive pocket magnifier and a flexible arm of some sort. The trouble with buying cheap magnifiers is that a lot of them are completely useless shapes. That is to say that they are neither spherical or hyperbolic lens profiles, but they're molded or ground to some freehand wavy biconic shape that makes them about as useful as a fishbowl. I found some glass pocket magnifiers on Ebay that have proven to be good. I just popped the rivet out of one and made a wire arm for it. It works well enough, but it tends to wobble if it's on a machine that vibrates. It needs a stiffer wire or some dampening ... or maybe I just need to replace the bad belts in the lathe. It's a promising candidate though.

Pretending to put a shear cut on some trashy porous cast aluminum bar

In the course of scouring china-mart for lens-shaped objects, I saw this thing and thought I had a brilliant idea. Increasing illumination helps a lot, partly because pupilary contraction increases depth of field. I could slap a USB power bank of some sort onto that and have a little portable magnifier light. Well, not really. The lenses in these things are absolutely useless. Luckily, it wasn't molded directly into the lamp.

Freehand wavy biconic lens-shaped object

I ended up taking another one of my glass pocket magnifiers and grinding the lens down to replace it. The gooseneck is barely able to hold the lamp up anyway, and the USB connector really is a lousy mechanical support. By the time I had the replacement lens fitted, I had pretty much condemned the USB power bank and the whole idea. I glued it to another magnet and figured I'll use it until I come up with something better.

It's easy to fix the lens. Just inherit an old lapidary wet grinder...

These aren't the best solutions, but they'll keep me busy for now. I plan on whipping up an articulated arm system in CAD. I'd like something simple that I can machine, but it might be nice to come up with a printable modular system.

If you're looking for a more substantial conclusion to this rambling post, there isn't one. These are not unique problems. Although it might be difficult for someone in good health to understand how frustrating and disempowering the experience of poor vision can be, plenty of people have it worse than I do.

Saturday, July 18, 2015

Animate a plot in MATLAB

As part of an ongoing exercise in pointlessly observing unavoidable disasters, I use a script to log the status of my internet connection via the information provided by the DSL modem. Normally, I just process this data into a pdf of daily plots of SNR and connection speed with average figures for each day. It's not too pretty, but it gets the point across easily enough.

For other aspects of the data, or maybe just to satisfy an idle whim, I thought it might be more interesting to visualize the data in motion. I had been playing with more complex figure manipulation recently, so I thought it would be a good challenge. I created a figure containing two axes which shows connection speed as well as signal and noise power. For simplicity, I used subplot_tight(), but there are other simplified methods of getting decent custom geometry out of figures with multiple axes. The background marking the days is generated simply:

% all the other data is prepared first
days=single(mod(round((t-12)/24),2)); % used to mark days
 
figure(1); clf; 
subplot_tight(4,1,2:4,0.05)
area(t,40*days-20,-20,'facecolor',0.9*[1 1 1],'edgecolor',[1 1 1]); hold on;
plot(t,PWRfill,'r:',t,NOISEfill,'b:',t,PWR,'r',t,NOISE,'b'); grid on; 
ax1=gca; 
set(gca,'layer','top'); % otherwise area() hides grid
xlabel('time (hours)'); ylabel('power (dBm)');
 
subplot_tight(4,1,1,0.05); 
area(t,1500*days,0,'facecolor',0.9*[1 1 1],'edgecolor',[1 1 1]); hold on;
plot(t,SPD); ax2=gca; set(gca,'layer','top'); grid on;
ylabel('Connection Speed kbps');
 
set(ax2,'XTickLabelMode', 'Manual','XTickLabel',[])
 
% keep subplots correlated on x when zooming
linkaxes([ax1 ax2],'x');

In the end, I decided against using the grid lines. They clutter the animation too much in my opinion. I left them in this code example simply to demonstrate the use of the 'Layer' property when using area() plots.

This also shows the overlay of duplicate datasets. The original data as logged contains thousands of gaps where the modem stopped reporting valid data, either because the connection was down or because the modem had locked up. While metrics like connection speed are only sensible to consider while the connection status is valid, things like noise power would really be most interesting at the extremes -- when the modem can't connect. I used inpaint_nans() to fill only the small gaps in the data, leaving longer gaps unfilled on the assumption that they represent modem lockups.

3-10-15 to 7-18-15 overall plot before animation

Once the figure is established with linked axes and everything is as desired, the animation routine incrementally steps through the data by setting the x-axis limits and capturing a frame as a .png file on disk. Matlab has an inbuilt function for generating movies (as an array of structs) and writing them as an .avi. In practice, this is a horrible burden on both ram and disk. For a video containing roughly 1300 frames, this method required a 4.9GB memory footprint and created an even larger uncompressed .avi file. The avifile() function can't utilize any compression in a linux environment, so it's pretty much useless to me. Caching images as png keeps memory use from scaling with frame count, and it allows external tools to do the video work. For some reason, Matlab does a much better job of creating png files than jpg files; they're much smaller too.

%% animate a plot
% assume axes are time-coordinated and linked
% acts generically only on one axis
 
width=36; % window width in hours
framestep=0.5; % time step between frames
outputdir='/home/assbutt/cad_and_projects/modemplots/';
 
% keep axes from autoscaling
a0=gca;
a=findall(gcf,'type','axes');
for n=1:1:length(a);
    axes(a(n));
    axis manual
end
axes(a0); % return to original axis
 
nframes=ceil((max(t)-width)/framestep)
perh=length(t)/max(t); % actual number of samples per hour
xlim=[0 width];
for f=1:1:nframes;
    set(gca,'xlim',xlim);
 
    % place a date label in the middle of the figure
    fdate=lineArray(round(perh*mean(xlim))); % date of center sample
    h=text(-100,0,sprintf('%s',fdate{:}),'horizontalalignment','center');
    set(h,'units','normalized','position',[0.5 -0.1 0],'fontsize',12);
 
    % image height is trimmed to an even pixel for video encoding
    pause(0.05); frame=getframe(gcf);
    frame.cdata=255-frame.cdata; sf=size(frame.cdata);
    imwrite(frame.cdata(1:end-(mod(sf(1),2)),:,:), ...
        [outputdir 'frame_' sprintf('%04d',f) '.png'],'png');
 
    xlim=xlim+framestep; delete(h);
end

I invert the image data during write because I normally use Matlab with the screen inverted. It's too easy for me to create plots that are difficult to read unless they stay that way.

The animation is a bit obnoxious because it grabs window focus constantly and makes it difficult to do anything else with the computer. I simply ran it on a second machine to avoid the hassle. Once everything was on disk, I compiled a video with avconv, though I suppose I could have gone directly to OpenShot with the image set.

avconv -r 24 -i frame_%04d.png -c:v libx264 -pix_fmt yuv420p out.mp4

I threw some extra titles and such on the ends using OpenShot and kicked it to YouTube so that it has somewhere to sit and collect dust.

Overall, I think the method works and helps get a point across. It's easy to see the transition in the noise characteristics over time. It's easy to identify the mid-day instability of the connection. I should have used a thicker linewidth when plotting things; scaled down, it gets difficult to see details and read text. Apologies for the blurriness.

Of course, maybe you don't want to lock axis scaling; maybe you want to change different aspects of the plot over the image sequence. The point is just to show that getframe() exists and avifile() is a fat ugly pig. Everything else is just hammering at graphics object properties in a loop.

I guess it's only appropriate

As I tried to explain last time, I consider this habit to be a unilateral affair. Until now, I'd no need to create a decoder. Out of curiosity, I decided to go ahead and make one; it was easy enough. To think all this time I hadn't known of base2dec(). I've been doing shit the hard way too many times. Oh well. Learning things is good.

%% text from unicode pictogram (arbitrary-base encoding)
clc; clear all; 
format compact
 
instring='▄▆▂▆ ▄█▄▄ ▄▆█▂ ▄▆█▂ ▄█▂█ ▄▆▆▂ ▄▆▆▄ ▄█▄▂';
charmap='▂▄▆█'; % length(charmap)=b
space=' '; % thin space
 
b=length(charmap);
for n=0:1:b-1;
    instring(findstr(instring,charmap(n+1)))=num2str(n);
end
 
instring(findstr(instring,space))=' ';
instring=str2num(instring); % simplifies indexing
outstring='';
 
for n=1:1:length(instring);
    outstring=[outstring char(base2dec(num2str(instring(n)),b))];
end
 
outstring

It works well enough. I suppose I could make it automagically find a minimal character set from the input pictogram string, but it wouldn't have any way of guessing the likely order. Granted, even a human would have to guess if the charmap used in encoding weren't conspicuously ordered.

Thursday, July 16, 2015

▒░▒▓▓ ▒░▒▓▒ ▒▒░▓░ ▒░▓▓▓

▄▂▂▄ ▄█▂█ ▂▆▂▂ ▄▆█▄ ▄█▄▄ ▄▆▂█ ▄▆▆▂ ▂▆▂▂ ▄▆▂▄ ▄█▂█ ▂▆▂▂ ▄▂▆▄ ▂▆▂▂ ▄▆█▄ ▄▆▂▄ ▄█▆▄ ▂▆▂▂ ▄▆▆▂ ▄▆▂▄ ▄█▄▆ ▄▆▄▄ ▂▆▂▂ ▄▆▄▄ ▄▆█▆ ▄▆▆▆ ▄▆██ ▄█▆▄ ▄▆▄▄ ▄▆▄▂ ▂▆▂▂ ▄▆▂█ ▄▆▄▄ ▄█▂▆ ▄█▄▂ ▄▆▂▄ ▄▆▆▄ ▄▆█▆ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▄█▂█ ▂▆▂▂ ▄▆▆▄ ▄▆█▆ ▂▆▂▂ ▄█▂▂ ▄▆▂▄ ▄█▂▆ ▄█▄▂ ▄▆▆▄ ▄▆▂█ ▄█▄▄ ▄▆█▂ ▄▆▂▄ ▄█▂▆ ▂▆█▂ ▂▆▂▂ ▄▂▆▄ ▂▆▄█ ▄█▄▆ ▄▆▄▄ ▂▆▂▂ ▄▆▆▂ ▄▆▂▄ ▄▆▄▂ ▂▆▂▂ ▄▆▂▄ ▂▆▂▂ ▄▆▄█ ▄█▂▆ ▄▆██ ▄█▄█ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆▂▂ ▄▆▄▂ ▄▆▆▄ ▄▆▄▆ ▄▆▄▆ ▄▆▆▄ ▄▆▂█ ▄█▄▄ ▄▆█▂ ▄█▄▂ ▄█▆▄ ▂▆▂▂ ▄▆▄▄ ▄▆█▆ ▄▆▆▆ ▄▆██ ▄█▆▄ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆▂▂ ▄▆▂▄ ▄▆█▆ ▄█▆▄ ▄█▄▂ ▄▆▆▂ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆▂▂ ▄▆▆▄ ▄▆█▆ ▂▆▂▂ ▄▆▄█ ▄▆▄▄ ▄▆█▆ ▄▆▄▄ ▄█▂▆ ▄▆▂▄ ▄▆█▂ ▂▆█▆ ▂▆▂▂ ▂▆▂▂ ▄▂▆▄ ▄█▄▂ ▂▆▂▂ ▄▆▂█ ▄█▂▆ ▄▆▄▄ ▄▆▄▄ ▄█▂▂ ▄█▂█ ▂▆▂▂ ▄▆▆▄ ▄▆█▆ ▄█▄▂ ▄▆██ ▂▆▂▂ ▄▆▄▄ ▄█▄▆ ▄▆▄▄ ▄█▂▆ ▄█▆▄ ▄█▄▂ ▄▆▆▂ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆█▂ ▂▆▂▂ ▄▆▄▆ ▄█▂▆ ▄▆██ ▄▆█▄ ▂▆▂▂ ▄▆█▄ ▄█▆▄ ▂▆▂▂ ▄▆▄▂ ▄█▂▆ ▄▆▄▄ ▄▆▂▄ ▄▆█▄ ▄█▂█ ▂▆▂▂ ▄█▄▂ ▄▆██ ▂▆▂▂ ▄▆█▄ ▄▆▄▄ ▄▆█▄ ▄▆██ ▄█▂▆ ▄█▆▄ ▂▆█▆ ▂▆▂▂ ▂▆▂▂ ▄▄▄█ ▄▆▆▄ ▄█▄▂ ▄▆▆▂ ▂▆▂▂ ▄█▄█ ▄▆▂▄ ▄▆█▆ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆▂▂ ▄▆▂█ ▄▆██ ▄▆█▆ ▄▆▄▆ ▄▆▆▄ ▄▆▄▂ ▄▆▄▄ ▄▆█▆ ▄▆▂█ ▄▆▄▄ ▂▆█▂ ▂▆▂▂ ▄▂▆▄ ▂▆▂▂ ▄▆▄▆ ▄▆▄▄ ▄▆▄▄ ▄▆█▂ ▂▆▂▂ ▄▆▂▄ ▄█▂█ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆██ ▄█▄▄ ▄▆▄█ ▄▆▆▂ ▂▆▂▂ ▄▂▆▄ ▂▆▂▂ ▄▆▆▂ ▄▆▂▄ ▄█▄▆ ▄▆▄▄ ▄▆█▆ ▂▆▄█ ▄█▄▂ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▄▄ ▂▆▂▂ ▄▆▂▆ ▄█▂▆ ▄▆▄▄ ▄▆▂▄ ▄█▄▂ ▄▆▆▂ ▂▆▂▂ ▄█▄▂ ▄▆██ ▂▆▂▂ ▄▆▄▄ ▄█▄▆ ▄▆▄▄ ▄▆█▆ ▂▆▂▂ ▄▆▄▄ ▄█▆▂ ▄█▂▂ ▄█▂▆ ▄▆▄▄ ▄█▂█ ▄█▂█ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▂▄ ▄█▄▂ ▂▆▂▂ ▄▂▆▄ ▂▆▂▂ ▄▆▂█ ▄▆▂▄ ▄▆█▆ ▂▆▄█ ▄█▄▂ ▂▆▂▂ ▄▆▂█ ▄▆▂▄ ▄█▂▆ ▄▆▄▄ ▂▆▂▂ ▄▆▂▆ ▄▆▄▄ ▄█▆▄ ▄▆██ ▄▆█▆ ▄▆▄▂ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▄▄ ▂▆▂▂ ▄█▂▆ ▄▆▄▄ ▄▆█▄ ▄▆▂▄ ▄▆▆▄ ▄▆█▆ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆▂▂ ▄▆▄▂ ▄▆▆▄ ▄█▂█ ▄▆▄█ ▄█▄▄ ▄█▂█ ▄█▄▂ ▂▆▂▂ ▄▆██ ▄█▂▆ ▂▆▂▂ ▄█▂▆ ▄▆▄▄ ▄▆▄█ ▄█▂▆ ▄▆▄▄ ▄█▄▂ ▂▆█▆

▄▄▂█ ▄█▄▂ ▄▆▆▄ ▄▆█▂ ▄▆█▂ ▂▆█▂ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▆▄ ▄█▂█ ▂▆▂▂ ▄▆▂▄ ▄▆█▆ ▄▆▆▄ ▄▆█▄ ▄▆▂▄ ▄▆█▂ ▂▆▂▂ ▄▆█▆ ▄▆▄▄ ▄▆▄▄ ▄▆▄▂ ▂▆▂▂ ▄█▄▂ ▄▆██ ▂▆▂▂ ▄▆▂█ ▄▆██ ▄▆█▄ ▄▆█▄ ▄█▄▄ ▄▆█▆ ▄▆▆▄ ▄▆▂█ ▄▆▂▄ ▄█▄▂ ▄▆▄▄ ▂▆▂▂ ▄█▄█ ▄▆▄▄ ▄▆█▂ ▄▆█▂ ▄█▂█ ▂▆▂▂ ▄▆▄▆ ▄█▂▆ ▄▆██ ▄▆█▄ ▂▆▂▂ ▄█▄█ ▄▆▆▄ ▄█▄▂ ▄▆▆▂ ▄▆▆▄ ▄▆█▆ ▂▆█▆ ▂▆▂▂ ▂▆▂▂ ▄▂▂▄ ▄█▂█ ▂▆▂▂ ▄▆▂▄ ▂▆▂▂ ▄▆█▄ ▄▆▂▄ ▄█▂█ ▄█▄▂ ▄█▄▄ ▄█▂▆ ▄▆▂▆ ▄▆▂▄ ▄█▄▂ ▄▆▆▄ ▄▆██ ▄▆█▆ ▂▆▂▂ ▄▆██ ▄▆▄▆ ▂▆▂▂ ▄▆█▆ ▄▆▄▄ ▄▆▂█ ▄▆▄▄ ▄█▂█ ▄█▂█ ▄▆▆▄ ▄█▄▂ ▄█▆▄ ▂▆█▂ ▂▆▂▂ ▄▂▆▄ ▂▆▂▂ ▄▆▂█ ▄▆██ ▄▆█▆ ▄█▄▂ ▄▆▆▄ ▄▆█▆ ▄█▄▄ ▄▆▄▄ ▂▆▂▂ ▄█▄▂ ▄▆██ ▂▆▂▂ ▄▆▄█ ▄▆██ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄█▂▆ ▄▆██ ▄█▄▄ ▄▆▄█ ▄▆▆▂ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▄▄ ▂▆▂▂ ▄▆█▄ ▄▆██ ▄█▄▂ ▄▆▆▄ ▄▆██ ▄▆█▆ ▄█▂█ ▂▆█▆ ▂▆▂▂ ▂▆▂▂ ▄▂▆▄ ▄▆▄▆ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▄▄ ▂▆▂▂ ▄▆▄▂ ▄█▂▆ ▄▆▆▄ ▄█▄▆ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆▂▂ ▄▆▂█ ▄▆██ ▄▆█▆ ▄▆▄▆ ▄▆█▂ ▄▆▆▄ ▄▆▂█ ▄█▄▂ ▂▆▂▂ ▄█▂▆ ▄▆▄▄ ▄▆█▄ ▄▆▂▄ ▄▆▆▄ ▄▆█▆ ▄█▂█ ▂▆▂▂ ▄▆▆▄ ▄▆█▆ ▄█▄▂ ▄▆▄▄ ▄█▂▆ ▄▆█▆ ▄▆▂▄ ▄▆█▂ ▂▆█▂ ▂▆▂▂ ▄▆▂▄ ▄▆█▆ ▄▆▄▂ ▂▆▂▂ ▄▆█▆ ▄▆██ ▂▆▂▂ ▄▆▄▄ ▄█▆▂ ▄█▄▂ ▄▆▄▄ ▄█▂▆ ▄▆█▆ ▄▆▂▄ ▄▆█▂ ▂▆▂▂ ▄▆▄▄ ▄█▆▂ ▄█▂▂ ▄▆▄▄ ▄▆▂█ ▄█▄▂ ▄▆▂▄ ▄█▄▂ ▄▆▆▄ ▄▆██ ▄▆█▆ ▄█▂█ ▂▆▂▂ ▄▆▄▄ ▄█▆▂ ▄▆▆▄ ▄█▂█ ▄█▄▂ ▂▆█▂ ▂▆▂▂ ▄▆▂▄ ▄▆█▂ ▄▆█▂ ▂▆▂▂ ▄█▂▂ ▄█▂▆ ▄▆▄▄ ▄█▄▂ ▄▆▄▄ ▄▆█▆ ▄█▂█ ▄▆▄▄ ▂▆▂▂ ▄▆▂▄ ▄█▂█ ▂▆▂▂ ▄▆▂█ ▄▆██ ▄▆█▄ ▄▆█▄ ▄█▄▄ ▄▆█▆ ▄▆▆▄ ▄▆▂█ ▄▆▂▄ ▄█▄▂ ▄▆▆▄ ▄▆██ ▄▆█▆ ▂▆▂▂ ▄▆█▄ ▄▆▂▄ ▄█▆▄ ▂▆▂▂ ▄▆▂▄ ▄█▂█ ▂▆▂▂ ▄█▄█ ▄▆▄▄ ▄▆█▂ ▄▆█▂ ▂▆▂▂ ▄▆▂▆ ▄▆▄▄ ▂▆▂▂ ▄▆▂▄ ▄▆▂▆ ▄▆▂▄ ▄▆█▆ ▄▆▄▂ ▄▆██ ▄▆█▆ ▄▆▄▄ ▄▆▄▂ ▂▆█▆

▄▂▂▄ ▄▆█▆ ▄▆▄▂ ▂▆▂▂ ▄█▂█ ▄▆██ ▂▆▂▂ ▄▆▆▄ ▄█▄▂ ▂▆▂▂ ▄▆▄█ ▄▆██ ▄▆▄▄ ▄█▂█ ▂▆▂▂ ▄▆▂▄ ▄█▄▂ ▂▆▂▂ ▄█▄▂ ▄▆▆▄ ▄▆█▄ ▄▆▄▄ ▄█▂█ ▂█▆█ ▂▆▂▂ ▂▆▂▂ ▄▆▂▄ ▄█▂█ ▂▆▂▂ ▄▆▄▄ ▄▆▂▄ ▄█▂█ ▄▆▆▄ ▄▆█▂ ▄█▆▄ ▂▆▂▂ ▄▆▂▄ ▄█▂█ ▂▆▂▂ ▄▆██ ▄▆█▆ ▄▆▄▄ ▂▆▂▂ ▄▆▂█ ▄▆▂▄ ▄▆█▆ ▂▆▂▂ ▄▆▆▂ ▄▆██ ▄▆█▂ ▄▆▄▂ ▂▆▂▂ ▄▆▄▂ ▄▆▆▄ ▄▆▂▄ ▄▆█▂ ▄▆██ ▄▆▄█ ▄█▄▄ ▄▆▄▄ ▂▆▂▂ ▄█▄█ ▄▆▆▄ ▄█▄▂ ▄▆▆▂ ▂▆▂▂ ▄▆▂▄ ▄▆█▆ ▂▆▂▂ ▄▆▄▄ ▄▆█▄ ▄█▂▂ ▄█▄▂ ▄█▆▄ ▂▆▂▂ ▄█▂▆ ▄▆██ ▄▆██ ▄▆█▄ ▂▆█▂ ▂▆▂▂ ▄▂▆▄ ▂▆▂▂ ▄▆▂▆ ▄█▂▆ ▄█▄▄ ▄█▂█ ▄▆▆▂ ▂▆▂▂ ▄▆▂▄ ▄█▂█ ▄▆▆▄ ▄▆▄▂ ▄▆▄▄ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▄▄ ▂▆▂▂ ▄▆▂▄ ▄█▂█ ▄█▂█ ▄█▄▄ ▄▆█▄ ▄█▂▂ ▄█▄▂ ▄▆▆▄ ▄▆██ ▄▆█▆ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▂▄ ▄█▄▂ ▂▆▂▂ ▄▆█▄ ▄█▆▄ ▂▆▂▂ ▄█▄█ ▄▆██ ▄█▂▆ ▄▆▄▂ ▄█▂█ ▂▆▂▂ ▄█▄█ ▄▆▆▄ ▄▆█▂ ▄▆█▂ ▂▆▂▂ ▄▆▂▆ ▄▆▄▄ ▂▆▂▂ ▄█▂▆ ▄▆▄▄ ▄▆▂▄ ▄▆▄▂ ▂▆█▆ ▂▆▂▂ ▂▆▂▂ ▄▂▂▄ ▄█▂█ ▂▆▂▂ ▄▆▂▄ ▄▆█▆ ▄█▆▄ ▂▆▂▂ ▄█▂▂ ▄█▂▆ ▄▆██ ▄█▂▂ ▄▆▂▄ ▄▆▄█ ▄▆▂▄ ▄▆█▆ ▄▆▄▂ ▄▆▆▄ ▄█▂█ ▄█▄▂ ▂▆▂▂ ▄█▄█ ▄▆▆▄ ▄▆█▂ ▄▆█▂ ▂▆▂▂ ▄▆▄▂ ▄▆▄▄ ▄▆█▂ ▄▆▆▄ ▄▆▄█ ▄▆▆▂ ▄█▄▂ ▂▆▂▂ ▄█▄▂ ▄▆██ ▂▆▂▂ ▄▆▄▂ ▄▆▄▄ ▄▆█▆ ▄█▆▄ ▂▆█▂ ▂▆▂▂ ▄▆█▄ ▄▆▂▄ ▄▆▆█ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆▂▂ ▄█▂█ ▄█▄▄ ▄█▂▆ ▄▆▄▄ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▂▄ ▄█▄▂ ▂▆▂▂ ▄█▄█ ▄▆██ ▄█▂▆ ▄▆▄▂ ▄█▂█ ▂▆▂▂ ▄▆▂▄ ▄█▂▆ ▄▆▄▄ ▂▆▂▂ ▄▆█▄ ▄▆▄▄ ▄▆▂▄ ▄▆█▆ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▄▆█▂ ▄▆▄▄ ▄█▂█ ▄█▂█ ▂▆▂▂ ▄▆▆▄ ▄█▂█ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▄▄ ▂▆▂▂ ▄▆▄▆ ▄▆▆▄ ▄█▂▆ ▄█▂█ ▄█▄▂ ▂▆▂▂ ▄█▂█ ▄█▄▂ ▄▆▄▄ ▄█▂▂ ▂▆▂▂ ▄█▄▂ ▄▆██ ▄█▄█ ▄▆▂▄ ▄█▂▆ ▄▆▄▂ ▂▆▂▂ ▄▆▄▄ ▄█▂█ ▄█▄▂ ▄▆▂▄ ▄▆▂▆ ▄▆█▂ ▄▆▆▄ ▄█▂█ ▄▆▆▂ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆▂▂ ▄▆██ ▄█▂▆ ▂▆▂▂ ▄▆▂▄ ▄▆▂█ ▄▆▆█ ▄▆█▆ ▄▆██ ▄█▄█ ▄▆█▂ ▄▆▄▄ ▄▆▄▂ ▄▆▄█ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▂▆▂▂ ▄█▄▂ ▄▆▆▂ ▄▆▄▄ ▂▆▂▂ ▄▆█▄ ▄▆▄▄ ▄▆▂▄ ▄▆█▆ ▄▆▆▄ ▄▆█▆ ▄▆▄█ ▄▆█▂ ▄▆▄▄ ▄█▂█ ▄█▂█ ▄▆█▆ ▄▆▄▄ ▄█▂█ ▄█▂█ ▂▆▂▂ ▄▆██ ▄▆▄▆ ▂▆▂▂ ▄▆▂▄ ▄▆▂█ ▄█▄▂ ▄▆▆▄ ▄▆██ ▄▆█▆ ▂▆█▆

% text to unicode pictogram (arbitrary-base encoding)
clc; clear all; 
format compact
 
instring=sprintf('And of course, I used Matlab to do it!');
charmap='▂▄▆█'; % length(charmap) determines numeric base
space=' '; % thin space
 
outstring='';
newchar='';
b=length(charmap);
bs=ceil(log(128)/log(b));
for n=1:1:length(instring);
    nc=double(instring(n));
    for m=1:1:bs;
        bitstring=floor(mod(nc,b^m)/b^(m-1));
        for k=0:1:bs-1;
            if bitstring==k
                newchar=[charmap(k+1) newchar];
            end
        end
    end
     
    outstring=[outstring newchar space];
    newchar='';
end
outstring

▄▂▂▄ ▄▆█▆ ▄▆▄▂ ▂▆▂▂ ▄▆██ ▄▆▄▆ ▂▆▂▂ ▄▆▂█ ▄▆██ ▄█▄▄ ▄█▂▆ ▄█▂█ ▄▆▄▄ ▂▆█▂ ▂▆▂▂ ▄▂▆▄ ▂▆▂▂ ▄█▄▄ ▄█▂█ ▄▆▄▄ ▄▆▄▂ ▂▆▂▂ ▄▂█▄ ▄▆▂▄ ▄█▄▂ ▄▆█▂ ▄▆▂▄ ▄▆▂▆ ▂▆▂▂ ▄█▄▂ ▄▆██ ▂▆▂▂ ▄▆▄▂ ▄▆██ ▂▆▂▂ ▄▆▆▄ ▄█▄▂ ▂▆▂▄

Saturday, July 11, 2015

Blend images with MATLAB

UPDATE: These scripts have been vastly improved. Find the current versions here.

In the course of doing some image mutilation in Matlab, I had need to perform some blend operations. Prior to that point, I had been using GIMP to perform the task, but I'd rather be able to automate things. There are probably a number of FEX submissions that could replace this function, but like usual, I decided to reinvent the wheel and write my own.

The following function makes available most of the common blend modes available in GIMP or Photoshop, but I decided to add a few extra things. Uncommon features include:

Adjustable dodge/burn amount (not just opacity!)
Several luminance-dependent analogs of common functions
Lightness and Intensity modes
Several hue permutation modes
Several color permutation (hue & saturation) modes
Inputs can be single images or image sequences (4-D inputs!)

The ability to handle 4-D inputs is a big plus for processing animation sequences. The function can take matching 3-D or 4-D array inputs, or a mismatched input of one image and a 4-D array. In the latter case, the single image is blended with all frames of the 4-D input. This makes it possible to easily blend a static overlay onto all frames of an animation, for instance.

This is probably in need of some work still; I'd like to add some sort of alpha support to things perhaps. In the meantime, it's working well. Many of the conversions use colorspace() from Pascal Getreuer.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

function  outpict=imblend(FG,BG,opacity,blendmode,amount)
%   IMBLEND(FG, BG, OPACITY, BLENDMODE,{AMOUNT})
%       blend images or imagesets as one would blend layers in GIMP or
%       Photoshop.  
%
%   FG, BG are RGB image arrays of same H,V dimension
%       both can be single images or 4-D imagesets of equal length
%       can also blend a single image with a 4-D imageset
%       mismatches of dimensions 3:4 result in array expansion
%       allows blending static overlays with an entire animation
%       mismatches of dimensions 1:2 are not supported
%   OPACITY is a scalar from 0 to 1
%       defines mixing of blended result and original BG
%   AMOUNT is a scalar (optional, default 1)
%       used to internally scale the influence of blend calculations
%   BLENDMODE is a string assignment (see list)
%
%
%   MODES: 
%       normal      
%       screen
%       overlay     (standard method)
%       softlight   (GIMP overlay)
%       hardlight
%       vividlight
%       hardmix     (similar to posterization)                  amount:[0 1]
%       posterize   (stronger influence from mask)
%       colordodge  (similar to GIMP dodge)                     amount:[0 1]
%       colorburn   (similar to GIMP burn)                      amount:[0 1]
%       lineardodge                                             amount:[0 1]
%       linearburn                                              amount:[0 1]
%       lighten RGB     (lighten only (RGB))
%       darken RGB      (darken only (RGB))
%       lighten Y       (lighten only (luma only))
%       darken Y        (darken only (luma only))
%       scale add       (add bg to fg deviation from mean)     amount:(-inf to +inf)
%       scale mult      (scale bg by mean-normalized fg)       amount:(-inf to +inf)
%       multiply
%       divide
%       addition
%       subtraction
%       difference
%       exclusion 
%       hue         
%       saturation  
%       value    
%       luma1       (uses colorspace() YIQ conversion)
%       luma2       (Image Processing toolbox YIQ conversion)
%       lightness   (approx identical to intensity)
%       intensity
%       color       
%       permute H>H     (rotates hue by mask hue)               amount:(-inf to +inf)
%       permute dH>H    (rotates hue by hue difference)         amount:(-inf to +inf)
%       permute Y>H     (rotates hue by mask luma)              amount:(-inf to +inf)
%       permute dY>H    (rotates hue by luma difference)        amount:(-inf to +inf)
%       permute H>HS     (rotates color by mask hue)            amount:(-inf to +inf)
%       permute dH>HS    (rotates color by hue difference)      amount:(-inf to +inf)
%       permute Y>HS     (rotates color by mask luma)           amount:(-inf to +inf)
%       permute dY>HS    (rotates color by luma difference)     amount:(-inf to +inf)
%
%   NOTE:
%       modes which accept 'amount' argument are marked with effective range
%       dH>H and dH>HS permutations are same as 'hue' when amount==-1
%       color permutations combine hue rotation and saturation blending
%           saturation blending is maximized when abs(amount)==1
%
%   CLASS SUPPORT:
%       Accepts images of 'uint8', 'double', and 'logical'
%       Return type is inherited from BG
%       In the case of a 'double' input, any image containing values >1
%       is assumed to have a white value of 255. 
 
%   SOURCES:
%       http://www.venture-ware.com/kevin/coding/lets-learn-math-photoshop-blend-modes/
%       http://www.deepskycolors.com/archive/2010/04/21/formulas-for-Photoshop-blending-modes.html
%       http://en.wikipedia.org/wiki/Blend_modes
%       http://en.wikipedia.org/wiki/YUV
%       http://www.kineticsystem.org/?q=node/13
%       http://www.simplefilter.de/en/basics/mixmods.html
 
if nargin ~= 5
    amount=1;
end
 
 
% i had intended to make this more class-insensitive, but i never need it
% output type is inherited from BG, assumes white value of either 1 or 255
inclassFG=class(FG);
inclassBG=class(BG);
if strcmp(inclassFG,'uint8')
    fgmax=255;
elseif strcmp(inclassFG,'double')
    if max(max(max(FG)))<=1
        fgmax=1;
    else
        fgmax=255;
    end
elseif strcmp(inclassFG,'logical')
    fgmax=1;
else
    disp('IMBLEND: unsupported class for FG')
    return
end
 
if strcmp(inclassBG,'uint8')
    bgmax=255;
elseif strcmp(inclassBG,'double')
    if max(max(max(BG)))<=1
        bgmax=1;
    else
        bgmax=255;
    end
elseif strcmp(inclassBG,'logical')
    bgmax=1;
else
    disp('IMBLEND: unsupported class for BG')
    return
end
 
 
% expand along dimension 3 where necessary
if size(FG,3)&ltsize(BG,3)
    FG=repmat(FG,[1 1 size(BG,3) 1]);
elseif size(FG,3)&gtsize(BG,3)
    BG=repmat(BG,[1 1 size(FG,3) 1]);
end
 
% check if height & width match
sFG=size(FG);
sBG=size(BG);  
if any(sFG(1:2)~=sBG(1:2)) 
    disp('IMBLEND: images of mismatched dimension')
    return
end
 
% check frame count and expand as necessary
if length(sFG)~=4 && length(sBG)~=4 % two single images
    images=1;
else
    if length(sFG)~=4 % single FG, multiple BG
        FG=repmat(FG,[1 1 1 sBG(4)]);
    elseif length(sBG)~=4 % multiple FG, single BG
        BG=repmat(BG,[1 1 1 sFG(4)]); sBG=size(BG);
    elseif sFG(4)~=sBG(4) % two unequal imagesets
        disp('IMBLEND: imagesets of unequal length')
        return
    end
    images=sBG(4);
end
 
% perform blend operations
outpict=zeros(sBG);    
for n=1:1:images
    I=double(BG(:,:,:,n))/bgmax;
    M=double(FG(:,:,:,n))/fgmax;
 
    switch lower(blendmode)
        case 'normal'
            R=M;
 
        case 'screen'
            R=1-((1-M).*(1-I));
 
        case 'overlay'  % actual standard overlay mode 
            hi=I>0.5; lo=I&lt=0.5;
            R=zeros(size(I));
            R(lo)=2*I(lo).*M(lo);
            R(hi)=1-2*(1-M(hi)).*(1-I(hi));
 
        case 'softlight' % same as GIMP 'overlay' due to legacy bug
            Rs=1-((1-M).*(1-I));
            R=(I.*((1-I).*M+Rs));
 
        case 'hardlight'
            hi=M&gt0.5; lo=M&lt=0.5;
            R=zeros(size(I));
            R(lo)=2*I(lo).*M(lo);
            R(hi)=1-2*(1-M(hi)).*(1-I(hi));
 
        case 'vividlight'  % test this; example formulae are inconsistent
            hi=M&gt0.5; lo=M&lt=0.5;
            R=zeros(size(I));
            R(lo)=1-(1-I(lo))./(2*M(lo));
            R(hi)=I(hi)./(1-2*(M(hi)-0.5));
 
        case 'posterize'  % actually a broken version of vividlight
            hi=M&gt0.5; lo=M&lt=0.5;
            R=zeros(size(I));
            R(lo)=(1-I(lo))./(2*(M(lo)-0.5));
            R(hi)=1-I(hi)./(1-2*M(hi));
 
        case 'hardmix' % ps mode similar to posterization
            amount=max(min(amount,1),0);
            Rs=M+I;
            R=Rs;
            R(Rs&gt1)=1*amount;
            R(Rs&lt1)=0;
 
        % DODGES/BURNS    
        case'colordodge'
            amount=max(min(amount,1),0);
            R=I./(1-M*amount);
 
        case 'colorburn'
            amount=max(min(amount,1),0);
            R=1-(1-I)./(M*amount+(1-amount));
 
        case 'lineardodge' % addition
            amount=max(min(amount,1),0);
            R=M*amount+I;
 
        case 'linearburn'
            amount=max(min(amount,1),0);
            R=M*amount+I-1*amount;
 
        % SIMPLE MATH OPS    
        case 'lighten rgb' % lighten only (RGB, no luminance)
            R=max(I,M);
 
        case 'darken rgb' % darken only (RGB, no luminance)
            R=min(I,M);
 
        case 'lighten y' % lighten only (based on luminance)
            Myiq=colorspace('RGB-&gtYIQ',M);
            Iyiq=colorspace('RGB-&gtYIQ',I);
            mask=Myiq(:,:,1)&gtIyiq(:,:,1);
            R=double(replacepixels(255*I,mask,255*M))/255;
             
        case 'darken y' % darken only (based on luminance)
            Myiq=colorspace('RGB-&gtYIQ',M);
            Iyiq=colorspace('RGB-&gtYIQ',I);
            mask=Myiq(:,:,1)&ltIyiq(:,:,1);
            R=double(replacepixels(255*I,mask,255*M))/255;
 
        case 'multiply'
            R=M.*I;
 
        case 'divide'
            R=I./(M+1E-3);
 
        case 'addition' % same as lineardodge
            R=M+I;
 
        case 'subtraction'
            R=I-M;
 
        case 'difference'
            R=abs(M-I);
 
        case 'exclusion'
            R=M+I-2*M.*I;
 
        case 'hue'
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,1)=Mhsv(:,:,1);
            R=hsv2rgb(Rhsv);    
 
        case 'saturation'
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,2)=Mhsv(:,:,2);
            R=hsv2rgb(Rhsv); 
 
        % V=max([R G B])
        % L=mean(max([R G B]),min([R G B]))
        % I=mean([R G B])
        % Y=[0.299 0.587 0.114]*[R G B]'
 
        case 'value'
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,3)=Mhsv(:,:,3);
            R=hsv2rgb(Rhsv); 
 
        % all colorspace() Y-swaps produce identical results within 1 LSB
        % (YUV, YIQ, YCbCr, YPbPr, YDbDr)
        case 'luma1' % swaps fg bg luma
            Myiq=colorspace('RGB-&gtYIQ',M);
            Ryiq=colorspace('RGB-&gtYIQ',I);
            Ryiq(:,:,1)=Myiq(:,:,1);
            R=colorspace('RGB&lt-YIQ',Ryiq);
 
        case 'luma2' % swaps fg bg luma (using IP toolbox)
            Myiq=rgb2ntsc(M);
            Ryiq=rgb2ntsc(I);
            Ryiq(:,:,1)=Myiq(:,:,1);
            R=ntsc2rgb(Ryiq); 
 
        % L and I swaps are calculated differently,  
        % but results are practically identical (within 1 LSB) 
        % for all available HSL and HSI conversion implementations
        case 'lightness' % swaps fg bg lightness
            Mhsl=colorspace('RGB-&gtHSL',M);
            Rhsl=colorspace('RGB-&gtHSL',I);
            Rhsl(:,:,3)=Mhsl(:,:,3);
            R=colorspace('RGB&lt-HSL',Rhsl);
 
        case 'intensity' % swaps fg bg intensity 
            Mhsi=colorspace('RGB-&gtHSI',M);
            Rhsi=colorspace('RGB-&gtHSI',I);
            Rhsi(:,:,3)=Mhsi(:,:,3);
            R=colorspace('RGB&lt-HSI',Rhsi);
 
        case 'color' % same as GIMP, swap H&S
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,1:2)=Mhsv(:,:,1:2);
            R=hsv2rgb(Rhsv); 
 
        % HUE PERMUTATIONS 
        case 'permute y>h' % permutes bg hue based on fg luma
            factors=[0.299 0.587 0.114];
            osize=size(M(:,:,1));
            cscale=repmat(reshape(factors,1,1,3),[osize 1]);
            Y=sum(M.*cscale,3);
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,1)=mod(Rhsv(:,:,1)+Y*amount,1);
            R=hsv2rgb(Rhsv); 
 
        case 'permute h>h' % permutes bg hue based on fg hue
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,1)=mod(Rhsv(:,:,1)+Mhsv(:,:,1)*amount,1);
            R=hsv2rgb(Rhsv);    
 
        case 'permute dy>h' % permutes bg hue based on luma difference
            factors=[0.299 0.587 0.114];
            osize=size(M(:,:,1));
            cscale=repmat(reshape(factors,1,1,3),[osize 1]);
            Ym=sum(M.*cscale,3);
            Yi=sum(I.*cscale,3);
            dY=Yi-Ym;
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,1)=mod(Rhsv(:,:,1)+dY*amount,1);
            R=hsv2rgb(Rhsv); 
 
        % note that dH+H permutation is same as a hue swap when amount==-1
        case 'permute dh>h' % permutes bg hue based on hue difference
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            dH=Rhsv(:,:,1)-Mhsv(:,:,1);
            Rhsv(:,:,1)=mod(Rhsv(:,:,1)+dH*amount,1);
            R=hsv2rgb(Rhsv); 
 
        % COLOR PERMUTATIONS (rotate hue and blend saturation)
        case 'permute y>hs' % permutes bg color based on fg luma
            factors=[0.299 0.587 0.114];
            osize=size(M(:,:,1));
            cscale=repmat(reshape(factors,1,1,3),[osize 1]);
            Y=sum(M.*cscale,3);
            amt=max(min(abs(amount),1),0); % needed since S-blending has limited range
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,1)=mod(Rhsv(:,:,1)+Y*amount,1);
            Rhsv(:,:,2)=amt*Mhsv(:,:,2)+(1-amt)*Rhsv(:,:,2);
            R=hsv2rgb(Rhsv); 
 
        case 'permute h>hs' % permutes bg color based on fg hue
            amt=max(min(abs(amount),1),0); % needed since S-blending has limited range
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,1)=mod(Rhsv(:,:,1)+Mhsv(:,:,1)*amount,1);
            Rhsv(:,:,2)=amt*Mhsv(:,:,2)+(1-amt)*Rhsv(:,:,2);
            R=hsv2rgb(Rhsv);    
 
        case 'permute dy>hs' % permutes bg color based on luma difference
            factors=[0.299 0.587 0.114];
            osize=size(M(:,:,1));
            cscale=repmat(reshape(factors,1,1,3),[osize 1]);
            Ym=sum(M.*cscale,3);
            Yi=sum(I.*cscale,3);
            dY=Yi-Ym;
            amt=max(min(abs(amount),1),0); % needed since S-blending has limited range
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            Rhsv(:,:,1)=mod(Rhsv(:,:,1)+dY*amount,1);
            Rhsv(:,:,2)=amt*Mhsv(:,:,2)+(1-amt)*Rhsv(:,:,2);
            R=hsv2rgb(Rhsv); 
 
        % note that dH+H permutation is same as a hue swap when amount==-1
        case 'permute dh>hs' % permutes bg color based on hue difference
            amt=max(min(abs(amount),1),0); % needed since S-blending has limited range
            Mhsv=rgb2hsv(M);
            Rhsv=rgb2hsv(I);
            dH=Rhsv(:,:,1)-Mhsv(:,:,1);
            Rhsv(:,:,1)=mod(Rhsv(:,:,1)+dH*amount,1);
            Rhsv(:,:,2)=amt*Mhsv(:,:,2)+(1-amt)*Rhsv(:,:,2);
            R=hsv2rgb(Rhsv); 
             
        % SCALE ADD treats FG as an additive gain map with a null point at its mean
        case 'scale add'
            Mstretch=imadjust(M,stretchlim(M));
            centercolor=mean(mean(Mstretch,1),2);
            R=zeros(size(I));
            for c=1:1:3;
                R(:,:,c)=I(:,:,c)+(Mstretch(:,:,c)-centercolor(:,:,c))*amount;
            end
 
        % SCALE MULT treats FG as a gain map with a null point at its mean
        case 'scale mult'
            Mstretch=imadjust(M,stretchlim(M));
            centercolor=mean(mean(Mstretch,1),2);
            R=zeros(size(I));
            for c=1:1:3;
                R(:,:,c)=I(:,:,c).*(Mstretch(:,:,c)./centercolor(:,:,c))*amount;
            end
             
        otherwise
            disp('IMBLEND: unknown blend mode');
            return
 
    end
 
    R=min(R,1); 
    R=max(R,0);
    R(isnan(R))=1;
    outpict(:,:,:,n)=bgmax*(opacity*R + I*(1-opacity));
end
 
outpict=cast(outpict,inclassBG);
 
return

I've been thinking about putting some polish on the rest of my image-garbling toolbox, so maybe I'll reveal some examples in time. I've also been toying with the idea of doing some animated plots of large datasets. That may be interesting.

Wednesday, July 1, 2015

Revisiting STARS-922 and Other Stuff

There seems to be some skepticism and curiosity that surround STARS-922 thermal compound (aka "heat sink plaster") and some similar products. I know I've wondered for a while myself and I've run across a handful of forum posts that question its purpose or properties. Places like DealExtreme and Ebay sell lots of this product along side other more common zinc-oxide thermal pastes with limited product differentiation. The poorly translated and dubious specifications and naming don't help to instill confidence that it is trustworthy of any particular application. On some occasions, customer reviews boast "strong" or "excellent conductivity" although i doubt either metric had been actually evaluated beyond the emotional satisfaction associated with the completion or failure of the parent project as determined by myriad other factors. I've found few satisfactorily detailed reviews, so I figured I'd add this bit of info to augment my halfass thermal testing from the prior post.

What the hell do you mean "plaster"?

Some of the questionable property specifications are:

Thermal conductivity: > 1.2W/m-K
Thermal Impedance: < 0.06 (what units?)
Clotting time: 3min (25 degree celsius) (all of it or just the exposed stuff?)
Strength of connected buildings: 25Kg (what?)
Temperature resistance: 200 degree celsius

Following some of my own experience and a few comments noting that it had little adhesive capability, I figured I would try to get a better feel for the mechanical properties of the product in a way that should at least lend intuition a better handle on what's appropriate to expect.

First off, this product is not a thermal grease. It is much more viscous and does indeed congeal on exposure to air. The cured product is slightly rubbery, but not as elastic as one might expect from a conventional RTV silicone caulking. The most common place I have seen this is in the bonding of the MCPCB to the radiator body of cheap LED retrofit lamps. In such an application, temperatures are low and mechanical stress is small. The MCPCB has only a small mass, and the aspect ratio ensures that it is difficult to produce a tension stress at the interface. This is all overshadowed by the advantage of being enclosed. Nothing can stress the bonded part except its own inertia. This seems to be a rather forgiving application, mechanically speaking. Will it work for larger SMT heat sinks or other things? We need to have a bit more understanding about its actual strength in shear and tension.

In order to test shear and tension strength, I put together a setup in the lathe consisting of various repurposed fixturing and a 1000# load cell with a selectable-gain instrumentation amplifier. Since I never finished the microcontroller DRO that was supposed to go with the amplifier I made for this load cell, I have to use the multimeter to take direct readings. Despite the extra clumsiness of operations, the results would be the same. To get the approximate peak value at break, I simply used the phone to record operations. Accepting that the periodic multimeter update is the result of an averaging operation and may lag the actual value, and forgiving my inability to produce a smooth and monotonically increasing tension by handscrew, the captured values should be acceptable for the intended purpose if averaged across a few samples.

some assembled samples and spring clamp for alignment

The shear samples consist of aluminum plates cut so that they produce a lap joint with approximately one square inch of area. The tension samples are old elevator bucket bolts that have had their head faces ground and scoured. All samples are ground and scoured in degreaser as needed to prepare a fresh surface each time. No particular assembly fixturing was used; no clamping was performed except to gently maintain alignment with adhesives that were flowable (e.g. molten hot glue). Samples were tested with heat sink plaster, some Permatex RTV silicone products, as well as generic superglue and two unknown hot glue products. One was a soft and flexible hobby-marketed product, and the other was a harder Arrow brand carpentry-marketed product, if I recall correctly.

In approximate terms, the heat sink plaster (HSP) is about on par with a low-strength RTV silicone. As with the prior thermal resistance testing, I ran a few sample groups of Permatex Ultra Grey RTV and High Temp Red RTV gasketing products. These represent the extremes of the product spectrum in terms of specified tensile strength and filler density. It's a fair bet to assume that if you'd be comfortable with the strength of your favorite RTV silicone, the strength of HSP will be sufficient for your task. There are a few things to consider, though.

One of the primary obstacles I ran into with this setup is that I'm testing air-cured materials in thin sections. With nonporous sample substrate, the geometry of the test patch does end up becoming important. Especially with attempts to cure products at room temperature, it was not uncommon to find samples that were only cured on an annulus around the perimeter of the bond patch. Most of these partially-cured sample groups were discarded for various reasons, typically nonuniformity or large variance within a group. In the case of some samples with uniform partial curing, I made my best attempt to measure the effective cured area using image processing after the break test. While I note that my goal here is to measure strength of the cured adhesive fraction, I did not make any differentiation in the thermal testing process. It may be worth knowing if the thermal resistance is much higher when the center of the sample patch has cured completely, but I can't exactly put the thermal test fixture in the oven with all its plastic bits. I'm also far too lazy to do larger test groups over a longer natural curing interval of say a month.

For what it's worth, the image processing used to calculate partial bond area is rather simple and goes to show that GIMP or Photoshop are useful for more than just ~~meandering artfaggotry~~ drawing and editing photos. It goes like this:

Remove uncured adhesive with solvent
Take photos of both halves
Transform image areas to geometric correspondence
Perform threshold operation to isolate cured material
Multiply images to get combined coverage
Find average pixel value

The example shown is not a particularly well-cured sample (42% cured), but it was from a group that were uniform.

In my attempts to expedite thorough curing of the samples, I tried several runs at elevated temperatures in the oven. Some of these resulted in reduced and typically inconsistent strength. In the case of the RTV, some of the extended cure times (96 hours at 70-80 °C) showed very low shear strength, and the bond patterns were not characteristic of radial curing. My guess is that as the silicone cured, it becomes increasingly difficult to outgas any remaining volatile products and pressure produced tear channels in the material near which further curing can take place.

some bad rtv shear samples

Some of the samples of heat sink plaster which were cured at high temperatures (100 °C) became dry and friable, taking on a tan color. These samples still exhibited normal shear strength, but had dramatically reduced tensile strength. Whether this loss of pliability and tensile strength occurs at lower temperatures over a longer time frame is unknown, but could potentially be a limiting factor in reliability. It also calls into direct question the product specifications. I doubt that the dry state is an intended condition for the product; consequently, I doubt that the 200 °C rating is realistic. The condition of samples cured at 70 °C was significantly better than those cured at only 30 °C more.

Overall, the product performance is more than adequate for small heat sinks, particularly low-height heat sinks. The thermal performance is not exceptional, but for the same area, it has a similar thermal resistance compared to a grease, despite forming a relatively thick interface film. I still would hesitate to refer to it as a self-shimming type TIM, as I doubt that it's intended to serve as an electrical insulator. That said, It's important to note the limitations of these tests.

All of the RTV silicone and HSP samples are very fresh and really don't represent thoroughly cured material. It may be reasonable to expect stronger bonds over a time span of a month or so, depending on the ability of the center of the bond patch to breathe. Similarly, these tests do not indicate whether strength or thermal performance will be sustained over much longer time frames at elevated temperatures. Longevity is still uncertain, and performance at temperatures beyond 100 °C is questionable.

That said, there are always other products available that fill a similar role. STARS-922 is simply the cheapest option. There are other cheap options out there in the realm of silicone products, and there are also epoxies for use when greater strength is desired and disassembly isn't intended. There are always the myriad products from Loctite and Dow for specific applications, though most of these are far too expensive for the hobbyist.

One other note as an aside: the hot glue samples were not applied with a glue gun. These, like most cases where I use hot glue, were applied with a heat gun. This is due to the simple fact that molten glue typically lacks the thermal capacity to heat a conductive substrate sufficiently to allow bonding to take place before the glue solidifies. Disregard for this simple fact appears widespread, as it seems every time I encounter hot glue in assembled products, it simply pops off the surfaces cleanly. So long as you can heat parts manually, hot glue does indeed have utility outside the realm of gluing googly eyes onto cotton balls. Even the cheap hobby products are incredibly strong so long as the application temperatures stay low.

Finally, I'll leave you with this bit of Matlab kludge that I put together to create the bar graph above. I always get so frustrated trying to use spreadsheet graphing utilities to do anyhing. It's impossible to get half the features that are desired, and the figures come out with cartoonish inelegance and poor repeatability. Granted, a configuration with proportionally-correspondent dual axes is a hassle even in Matlab, but at least Matlab provides a superabundance of flexibility. One might ask why anyone would ever need dual x or y axes, but the simple answer seems like a terribly common necessity to me: depicting one set of data in both metric and imperial units.

%% plot adhesive properties
clc; clf; clear all
format compact;
 
xt={'HT RTV 48h @ 70dC' 'HSP 24h @ 25dC' 'HSP 48h @ 70dC' 'UG RTV 48h @ 70dC' ...
    'Super Glue (generic)' 'Hot Glue (soft clear)' 'Hot Glue (hard amber)'};
ytkpa=[1299 1366 1860 2418 2552 3354 4758];
scalefactor=690/4758;
ytpsi=ytkpa*scalefactor;
scalelabels={'kPa' 'PSI'};
titlestring='Tensile Stress at Break';
ystretch=0; % stretches bar widths for better fit
 
% create first axis
hl1 = barh(1:1:length(xt),ytkpa);
ax1 = gca;
set(ax1,'XColor','r','YColor','k','box','off','color','none','YTickLabel',xt,...
    'XAxisLocation','bottom','YAxisLocation','left','XGrid','on',...
    'TickLength',[0.015 0],'ylim',[0+ystretch length(xt)+1-ystretch])
 
% this hackery removes the top and right tick marks from ax1
b = axes('Position',get(ax1,'Position'),'box','on','xtick',[],'ytick',[]);
axes(ax1); linkaxes([ax1 b]);
 
% create second axis
ax2 = axes('Position',get(ax1,'Position'),'XColor','b','YColor','k',...
    'XAxisLocation','top','YAxisLocation','right','Color','none',...
    'XGrid','on','Ytick',[],'TickLength',[0.02 0],...
    'xlim',get(ax1,'xlim')*scalefactor,'ylim',get(ax1,'ylim'));
hold(ax2,'on')
hl2 = barh(1:1:length(xt),ytpsi,'Parent',ax2);
set(hl2,'facecolor',[1 1 1]*0.5,'edgecolor',[1 1 1]*0);
set(get(hl2,'BaseLine'),'color','k')
 
% set and position the labels and titles
xl1=xlabel(ax1,scalelabels(1),'units','normalized','position',[-0.1 -0.02 0]);
xl2=xlabel(ax2,scalelabels(2),'units','normalized','position',[-0.1 1 0]);
t1=title(ax2,titlestring,'units','normalized','position',[0.5 1.09 0],'fontweight','bold');
 
% apply value labels on bars (not exactly robust geometry retention)
for n=1:1:length(xt)
    scalemax=max(get(ax2,'xlim'));    
    valb=text(ytpsi(n)-0.01*scalemax,n-0.0,['\color {red}' num2str(round(ytkpa(n))) ...
        ' \color {blue}' num2str(round(ytpsi(n)))], ...
        'horizontalalignment','right','fontweight','bold');
end
 
% this mess massages the plot to fit well within its container
marginoffset=[0.04 0.01];
margins=get(ax2,'tightinset'); h=[margins(1) margins(3)]; v=[margins(2) margins(4)];
set(ax1,'units','normalized','outerposition',[[h(1) v(1)]+marginoffset 1-sum(h) 1-sum(v)])
set(ax2,'units','normalized','outerposition',[[h(1) v(1)]+marginoffset 1-sum(h) 1-sum(v)])
set(b,'units','normalized','outerposition',[[h(1) v(1)]+marginoffset 1-sum(h) 1-sum(v)])
 
% updating axis sync without this is a bunch of shit
% similar to 'linkaxes()', but works with unequal axis limits (and 3 axes)
addlistener(ax1,'Position','PostSet',@(src,evt) updateAxis(ax1,ax2,b));

function updateAxes(a1,a2,a3)
    % callback function for updating internal plot geometry
    % to keep axis objects aligned if container geometry is altered
    xLim1=get(a2,'xLim');
    yLim1=get(a2,'yLim');
 
    set(a2,'position',get(a1,'position'));
    set(a3,'position',get(a1,'position'));
    set(a2,'xLim',xLim1);
    set(a2,'yLim',yLim1);
end