function out = cimage(cim, lim, cs)
%CIMAGE displays a complex valued data set as RGB, representing phase by
% color and amplitude by luminance.
%
% cimage(cim) displays the complex valued 2D matrix cim using autoscaling.
%
%
% cimage(cim, lim) displays the complex valued 2D matrix cim using the
%   two-element vector lim as lower and upper limits of the amplitude.
%
% cimage(cim, lim, cs) uses the color space cs, either 'hsv' (default),
%   'hsvwhite' or 'husl'
%
% out = cimage(...) outputs the RGB image to out instead of displaying it
%
%
%Example (creates a color bar/circle for reference):
% [re, im] = meshgrid(linspace(-1,1,1000));
% c = re + 1i*im;
% c(abs(c)>1) = NaN;
% cimage(c, [0 1])
%
%See also HSV2RGB, HUSL2RGB

%Revision History:
%120301: Created (Benjamin Judkewitz)
%120921: Updated help (Benjamin Judkewitz)
%160217: Enabled CLim changes and updates


    if nargin < 3 || isempty(cs)
        cs = 'hsv';
    end

    %if no limits are specified, normallize to maximum amplitude
    if nargin < 2 || isempty(lim)
        lim = [0 max(abs(cim(isfinite(cim(:)))))];
        disp(['lim: ' mat2str(lim)]);
    end

    %if only one limit is specified, assume it's the upper limit
    if numel(lim) == 1
        lim = [0 lim];
    end


    %this does the main work
    %rgb = complex2rgb(cim, lim);

    %display image if no output argument is provided - otherwise output data
    if nargout == 0
        hIm = image(nan([size(cim) 3]));
        hAx = gca;
        lh = addlistener(hAx,'CLim','PostSet',@(varargin) set(hIm, 'CData', complex2rgb(cim, hAx.CLim, cs))); %hAx.UserData = lh;
        hIm.DeleteFcn = @(varargin) delete(lh);
        hAx.CLim = lim;
    else
        out = complex2rgb(cim, lim, cs);
    end


function rgb = complex2rgb(cim, lim, method)
    switch (method)
        case 'hsv'
            %prepare H,S,V components
            H = angle(cim)/2/pi+.5;
            S = ones(size(cim));
            V = abs(cim);

            %scale V according to lim:
            V = (V-lim(1))/diff(lim);
            %S(V>1) = 0.2; %make them near-white
            V(V>1) = 1;
            V(V<0) = 0;

            %for white as background uncommment the following line to swap S and V:
            %temp = V; V = S; S = temp;

            %concatenate HSV and convert to RGB
            chsv = cat(3, H,S,V);
            rgb = zeros(size(chsv));
            for i = 1:size(chsv,4)
                rgb(:,:,:,i) = hsv2rgb(chsv(:,:,:,i));
            end

        case 'hsvwhite'
            H = angle(cim)/2/pi+.5;
            S = ones(size(cim));
            V = abs(cim);

            %scale V according to lim:
            V = (V-lim(1))/diff(lim);
            %S(V>1) = 0.2; %make them near-white
            V(V>1) = 1;
            V(V<0) = 0;

            temp = V; V = S; S = temp;

            %concatenate HSV and convert to RGB
            chsv = cat(3, H,S,V);
            rgb = zeros(size(chsv));
            for i = 1:size(chsv,4)
                rgb(:,:,:,i) = hsv2rgb(chsv(:,:,:,i));
            end
        case 'husl'
            %husl = (husl2rgb(shiftdim([linspace(0,360, 256)' ones(256, 1)*100 ones(256, 1)*60], -1))); mat2str(shiftdim(husl, 1))
            husl = [1 0.3004 0.5048;1 0.3035 0.4884;1 0.3065 0.4713;1 0.3095 0.4533;1 0.3125 0.4343;1 0.3154 0.4141;1 0.3184 0.3924;1 0.3213 0.369;1 0.3242 0.3434;1 0.3272 0.315;1 0.3301 0.2827;1 0.3331 0.2449;1 0.336 0.1981;1 0.339 0.1331;0.9998 0.3425 0;0.9814 0.3629 0;0.9644 0.3803 0;0.9484 0.3953 0;0.9335 0.4086 0;0.9194 0.4203 0;0.9061 0.4308 0;0.8935 0.4404 0;0.8816 0.449 0;0.8701 0.4569 0;0.8592 0.4642 0;0.8488 0.471 0;0.8388 0.4772 0;0.8291 0.4831 0;0.8198 0.4885 0;0.8108 0.4937 0;0.8021 0.4985 0;0.7936 0.5031 0;0.7854 0.5074 0;0.7774 0.5116 0;0.7696 0.5155 0;0.7619 0.5193 0;0.7544 0.5229 0;0.7471 0.5263 0;0.7399 0.5297 0;0.7328 0.5329 0;0.7258 0.536 0;0.7189 0.539 0;0.7121 0.5419 0;0.7054 0.5447 0;0.6987 0.5474 0;0.6921 0.5501 0;0.6855 0.5527 0;0.679 0.5552 0;0.6724 0.5577 0;0.6659 0.5602 0;0.6594 0.5626 0;0.6529 0.5649 0;0.6464 0.5672 0;0.6398 0.5695 0.0001294;0.6332 0.5718 0.0004896;0.6266 0.574 0.000862;0.6199 0.5762 0.001247;0.6132 0.5784 0.001645;0.6063 0.5805 0.002056;0.5994 0.5827 0.002481;0.5924 0.5848 0.002921;0.5853 0.5869 0.003375;0.5781 0.5891 0.003845;0.5707 0.5912 0.00433;0.5632 0.5933 1.479e-05;0.5555 0.5954 0;0.5476 0.5976 0;0.5395 0.5997 0;0.5311 0.6018 0;0.5225 0.604 0;0.5136 0.6062 0;0.5044 0.6084 0;0.4948 0.6106 0;0.4848 0.6128 0;0.4744 0.6151 0;0.4635 0.6174 0;0.452 0.6197 0;0.4399 0.6221 0;0.4271 0.6245 0;0.4134 0.627 0;0.3988 0.6295 0;0.383 0.632 0;0.3659 0.6347 0;0.3472 0.6374 0;0.3266 0.6401 0;0.3034 0.643 0;0.2769 0.6459 0;0.2459 0.6489 0;0.2079 0.6521 0;0.1577 0.6553 0;0.07519 0.6587 0;0.01933 0.659 0.09655;0.01814 0.658 0.1674;0.01703 0.6569 0.2133;0.01599 0.656 0.2483;0.01502 0.655 0.277;0.01411 0.6542 0.3014;0.01325 0.6533 0.3227;0.01245 0.6525 0.3416;0.0117 0.6517 0.3586;0.01099 0.651 0.374;0.01033 0.6503 0.3882;0.009701 0.6496 0.4012;0.009112 0.6489 0.4134;0.008558 0.6483 0.4247;0.008037 0.6476 0.4353;0.007546 0.647 0.4454;0.007085 0.6464 0.4549;0.006652 0.6459 0.4639;0.006246 0.6453 0.4724;0.005865 0.6447 0.4806;0.005508 0.6442 0.4885;0.005175 0.6437 0.496;0.004865 0.6431 0.5033;0.004576 0.6426 0.5103;0.004309 0.6421 0.517;0.004062 0.6416 0.5236;0.003835 0.6411 0.5299;0.003628 0.6406 0.5361;0.003439 0.6401 0.5422;0.00327 0.6397 0.5481;0.003118 0.6392 0.5538;0.002985 0.6387 0.5595;0.00287 0.6382 0.565;0.002772 0.6378 0.5705;0.002692 0.6373 0.5758;0.002629 0.6368 0.5811;0.002584 0.6363 0.5863;0.002556 0.6359 0.5914;0.002545 0.6354 0.5965;0.002552 0.6349 0.6016;0.002576 0.6345 0.6066;0.002618 0.634 0.6116;0.002678 0.6335 0.6165;0.002756 0.633 0.6215;0.002853 0.6325 0.6264;0.002968 0.632 0.6314;0.003103 0.6315 0.6363;0.003257 0.631 0.6412;0.003432 0.6305 0.6462;0.003627 0.63 0.6512;0.003843 0.6294 0.6562;0.004082 0.6289 0.6612;0.004343 0.6283 0.6663;0.004629 0.6278 0.6715;0.004938 0.6272 0.6767;0.005274 0.6266 0.6819;0.005636 0.626 0.6873;0.006027 0.6254 0.6927;0.006448 0.6247 0.6982;0.0069 0.624 0.7038;0.007385 0.6233 0.7096;0.007905 0.6226 0.7154;0.008462 0.6219 0.7214;0.009059 0.6211 0.7276;0.009698 0.6203 0.7339;0.01038 0.6195 0.7404;0.01112 0.6186 0.7471;0.0119 0.6177 0.754;0.01275 0.6168 0.7612;0.01365 0.6157 0.7686;0.01462 0.6147 0.7763;0.01567 0.6136 0.7843;0.0168 0.6124 0.7926;0.01801 0.6111 0.8014;0.01933 0.6097 0.8106;0.02075 0.6083 0.8202;0.0223 0.6067 0.8304;0.02398 0.6051 0.8411;0.02582 0.6032 0.8525;0.02784 0.6013 0.8646;0.03005 0.5991 0.8776;0.0325 0.5967 0.8915;0.03522 0.5941 0.9065;0.03826 0.5911 0.9227;0.04163 0.5878 0.9404;0.04526 0.5841 0.9597;0.04918 0.5798 0.981;0.08113 0.5748 1;0.1867 0.5689 1;0.2495 0.563 1;0.2976 0.5571 1;0.3376 0.5511 1;0.3723 0.5451 1;0.4034 0.539 1;0.4316 0.5329 1;0.4577 0.5267 1;0.4819 0.5204 1;0.5047 0.514 1;0.5263 0.5075 1;0.5469 0.5009 1;0.5666 0.4942 1;0.5855 0.4872 1;0.6039 0.4801 1;0.6216 0.4728 1;0.6389 0.4653 1;0.6558 0.4576 1;0.6723 0.4496 1;0.6885 0.4413 1;0.7044 0.4327 1;0.7201 0.4238 1;0.7357 0.4145 1;0.751 0.4047 1;0.7663 0.3944 1;0.7815 0.3836 1;0.7966 0.3721 1;0.8117 0.3599 1;0.8269 0.3468 1;0.842 0.3327 1;0.8572 0.3175 1;0.8726 0.3008 1;0.888 0.2823 1;0.9036 0.2616 1;0.9194 0.238 1;0.9354 0.2102 1;0.9517 0.1763 1;0.9682 0.1313 1;0.9851 0.05601 1;0.9938 0.008448 0.9933;0.9981 0.007782 0.9807;0.9996 0.03717 0.9658;0.9997 0.07216 0.95;0.9998 0.09601 0.9348;0.9999 0.1146 0.92;0.9999 0.13 0.9056;1 0.1432 0.8916;1 0.1548 0.8779;1 0.1651 0.8644;1 0.1744 0.8513;1 0.1829 0.8384;1 0.1907 0.8256;1 0.198 0.8131;1 0.2047 0.8007;1 0.211 0.7884;1 0.217 0.7763;1 0.2226 0.7642;1 0.228 0.7522;1 0.233 0.7402;1 0.2379 0.7283;1 0.2425 0.7164;1 0.247 0.7044;1 0.2513 0.6924;1 0.2554 0.6803;1 0.2594 0.6682;1 0.2633 0.6559;1 0.2671 0.6435;1 0.2707 0.631;1 0.2743 0.6182;1 0.2778 0.6052;1 0.2812 0.592;1 0.2845 0.5785;1 0.2878 0.5646;1 0.291 0.5504;1 0.2942 0.5357;1 0.2973 0.5205;1 0.3004 0.5048];
            rgb = interp1(linspace(-pi, pi, 256)', husl, angle(cim));
            a = (abs(cim)-lim(1)) / diff(lim);
            a(a>1) = 1;
            a = repmat(a, [1 1 3]);
            rgb = a .* rgb;
    end
