% [op 24th June 21] algorithm to recolour images; the aim here is to change the average chromaticity difference % of the images; please see Supplementary Material, page 9, of the paper % Visual discomfort and variations in chromaticity in art and nature for a % description of the algorithm % ========================================= % colour transforms etc. % ========================================= % ===== get colour toolbox ===== addpath('C:\Users\op5\Documents\MATLAB3\other_toolboxes3\CompColorScience') % ===== prepare transforms ===== cform_srgb2xyz = makecform('srgb2xyz'); % from srgb to CIE XYZ cform_xyz2upvpl = makecform('xyz2upvpl'); % from CIEXYZ to CIE LUV (aka CIE 1976 L*, u*, v*, and ??? CIE UCS 1976 as in Haigh Wilkins VR 2013) cform_xyz2srgb = makecform('xyz2srgb'); % to srgb from CIE XYZ cform_upvpl2xyz = makecform('upvpl2xyz'); % to CIEXYZ from CIE LUV % ========================================== % spectrum locus in CIELUV % ========================================== % ===== load basic spectrum locus ===== load('basic_upvp_spectrum_locus_v2.mat', 'bordupvp') % figure, plot(bordupvp(:,1)', bordupvp(:, 2)') xSL = bordupvp(:,1)'; ySL = bordupvp(:,2)'; % ===== get centre of space ===== (will be the centre of all rotations below) uvC = applycform([1 1 1]/3, cform_xyz2upvpl); % ===== we need to know the distance between the spectrum locus and the % centre for a dense set of directions ===== Ndir = 36; % re-center to uvC xSL = xSL - uvC(1); ySL = ySL - uvC(2); % figure, plot(xSL, ySL) % convert to polar coordinates [thetaSc, rhoSc] = cart2pol(xSL, ySL); % rhoSc is the scaling factor for thetaSc [thetaScs, Ord] = sort(thetaSc); rhoScs = rhoSc(Ord); % remove the duplicated points for interp1 (NB: points were duplicated to have full outline graphically) % [C,ia,ic] = unique(thetaScs); % look at ic % length(thetaScs) - length(unique(thetaScs)) thetaScs(19) = []; rhoScs(19) = []; figure, polarplot(thetaScs, rhoScs) % to make sure we don't miss anything we add a point at pi and -pi thetaScs = [-pi, thetaScs, pi]; rhoEnd = 0.5*(rhoScs(1) + rhoScs(end)); rhoScs = [rhoEnd, rhoScs, rhoEnd]; figure, polarplot(thetaScs, rhoScs) % ======================================== % image set % ======================================== % ===== choose a source folder and type of image (extension) ===== sourceFolder = 'choose source folder (flder continaing the images to be transformed)'; ext = '*.png'; % ===== get corresponding structure ===== imList = dir([sourceFolder, ext]); nbIm = size(imList, 1); % ========================================== % define set of transformations % ========================================== % simple rotations Phi1s = [-pi/2, -pi/4, 0, pi/4, pi/2]; % first rotation Phi2s = [-pi/2, -pi/4, 0, pi/4, pi/2]; % second rotation % families families = {'power','pseudo-logistic'}; % exponents lambdas = [0.01, 0.5, 4, 7, 12, 16, 20]; % N.B.: complexity is 5x5x2x7 = 350 % =========================================== % image destination % =========================================== destFolder = 'choose destination folder'; % =============================================== % loop over images > loop over transforms % =============================================== total = 0; for ii = 1:nbIm % read image imName = imList(ii).name; im = imread([sourceFolder, imName]); % save as reference image (i.e., the starting, or untransformed image) imwrite(im, [destFolder imName(1:end- 4) '_ref.png']) % to upvp im_xyz = applycform(im, cform_srgb2xyz); im_upvpl = applycform(im_xyz, cform_xyz2upvpl); % NB: l is u = reshape( im_upvpl(:,:,1), 1, []); v = reshape( im_upvpl(:,:,2), 1, []); % figure, plot(u, v, 'ok') % re-center to uvC/send to polar u = u - uvC(1); v = v - uvC(2); [thetaIm, rhoIm] = cart2pol(u, v); % figure, polarplot(thetaScs, rhoScs), hold on, polarplot(thetaIm, rhoIm, 'ok') for jj = 1:length(Phi1s) % loop over first rotations Phi1 = Phi1s(jj); thetaIm2 = thetaIm + Phi1; thetaIm2(thetaIm2 > pi) = thetaIm2(thetaIm2 > pi) - 2*pi; thetaIm2(thetaIm2 < -pi) = thetaIm2(thetaIm2 < -pi) + 2*pi; for kk = 1:length(families) % loop over families family = families{kk}; for ll = 1:length(lambdas) % loop over parameters within families lambda = lambdas(ll); % strech and compress (we use a monotonic remapping of [-pi, pi]) thetaIm3 = monotonic_remapping_interval(thetaIm2, -pi, pi, family, lambda); for mm = 1:length(Phi2s) % loop applying the composed transformations Phi2 = Phi2s(mm); % display infor for debugging and information str = 'current parameters are: image id-%d, angle1-%f, family-%s, exponent-%f, angle2-%f.\n'; fprintf(str, ii, 180*Phi1/pi, family, lambda, 180*Phi2/pi); % keep distance ratio to spectrum locus thetaIm4 = thetaIm3 + Phi2; thetaIm4(thetaIm4 > pi) = thetaIm4(thetaIm4 > pi) - 2*pi; thetaIm4(thetaIm4 < -pi) = thetaIm4(thetaIm4 < -pi) + 2*pi; rad1 = interp1(thetaScs, rhoScs, thetaIm); rad2 = interp1(thetaScs, rhoScs, thetaIm4); kappa = rad2./rad1; rhoIm4 = rhoIm.*kappa; % plot output figure, polarplot(thetaScs, rhoScs), hold on, polarplot(thetaIm, rhoIm, 'ok'), hold on, polarplot(thetaIm4, rhoIm4, 'or') % send to cartesian/de-center to uvC [u4, v4] = pol2cart(thetaIm4, rhoIm4); u4 = u4 + uvC(1); v4 = v4 + uvC(2); u4 = reshape(u4, size( squeeze(im_upvpl(:,:,3)))); v4 = reshape(v4, size( squeeze(im_upvpl(:,:,3)))); im2upvpl = zeros(size(im_upvpl)); im2upvpl(:,:, 1)= u4; im2upvpl(:,:, 2)= v4; im2upvpl(:,:, 3)= im_upvpl(:,:,3); total = total + 1; sum(kappa(isinf(kappa))) if total == 6 % just for debugging/testing stop_here = 1; end im2xyz = applycform(im2upvpl, cform_upvpl2xyz); % NB: l is the third dimension! im2 = applycform(im2xyz, cform_xyz2srgb); % % show output if needed % tempFig = figure; % subplot 121, imshow(im) % subplot 122, imshow(im2) % pause(1.1) % close all % save output imwrite(im2, [destFolder imName(1:end- 4) ' _angle1_' num2str(180*Phi1/pi) '_family_' family '_exponent_' num2str(lambda) '_angle2_' num2str(180*Phi2/pi) '.png']) end end end end end