% calibrate.m % % Depository for random code, until I figure out how to partition % it better. % % These values depend on the set of images % TakeTwo/cam1 saturationValue = 180; minValue = 10; macBethSpectrum='TakeTwo/cam1/t2mbSpecIndiv.mat'; imageDir='TakeTwo/cam1'; % Load macBeth color checker photospectrometer data % This loads mbSpec, the intensity per wavelength for each % checker, and mbWave, the wavelengths at which the intensities % were measured. The checkers are ordered by row, starting with % the grayscale row first. Here are a few for reference: % Checker 1: white 6: black % checker 7: blue 8: green 9: red load(macBethSpectrum); load XYZ; % This XYZ has values for 370:730 % Spectrometer output is 380:780 % Grab just the portion of XYZ that we need (380-728) % and portion of spectrometer output to match XYZ input range neededIndices = [11:4:359]; cropXYZ = XYZ(neededIndices,:); macSpec = mbSpec(1:88,:); macWave = mbWave(1:88); % Determine XYZ values of color checkers, % then scale Y values so white square has Y of 1.0 macXYZ = cropXYZ' * macSpec; macXYZ = macXYZ / macXYZ(2,1); % plot(mbWave,mbSpec(:,1)) % plot the white spectra -- reddish light! % plot(mbWave,mbSpec(:,7)) % plot the blue spectra % Test this against the photospectrometer output (correct). % macChroma = chromaticity(macXYZ); % Load the images. % Grabbing the color patches out of each image set is a bit % different, so I've put a function called "grabEachPatch" in the % directories of each set of images. It averages the dark images, % subtracts that value from the macBeth images, cd(imageDir); measRGB=grabEachPatch; cd ../.. % Doh! This exercise shows two things--first, that the sensor is % saturating, and second, it's response is linear until it saturates. % Redo measurements with non-saturating values. % This is a function of the lighting, too--the red sensor was % pegging first on the grayscale color checkers. figure; hold on plot(macXYZ(2,1:6),measRGB(1,1:6),'ro-'); plot(macXYZ(2,1:6),measRGB(2,1:6),'go-'); plot(macXYZ(2,1:6),measRGB(3,1:6),'bo-'); title('Measured R, G and B vs. measured Y for gray scale checkers'); xlabel('Luminance of checker;); ylabel('Pixel intensity'); hold off % Pick out the color checkers that aren't saturated or too low. maxVals = max(measRGB); minVals = min(measRGB); validCheckers = find(maxVals < saturationValue & minVals > minValue) validXYZ = macXYZ(:,validCheckers) validRGB = measRGB(:,validCheckers) %%% MODEL FITTING %%% Computing errors to data used for fit. % Ok, now it's finally time to do some model fitting. Start with % a 3x3 transformation from RGB to XYZ. See how good the fit % is. Does it matter which direction we go with these % transformations? I'm more interested in going from RGB back to % XYZ. % First Order Model, 3x3 Transform % validRGB'*colorTrans' = validXYZ' colorTrans = validRGB'\validXYZ'; colorTrans = colorTrans' sensorXYZ=colorTrans*validRGB; mse1=3*mse(sensorXYZ,validXYZ) %This is mean squared Euclidean distance in XYZ space. % Now let's calculate the error is CIE-Lab space. % Use the measured white point from the photospectrometer. whitePoint = macXYZ(:,1); validLab=xyz2lab(validXYZ',whitePoint'); measLab=xyz2lab(sensorXYZ',whitePoint'); mseLab1=3*mse(measLab,validLab) % Second fitting method % R G B RG GB BR [colorTrans2] = Model2(validRGB,validXYZ); [rgb2] = Model2Rgb(validRGB); sensorXYZ2=colorTrans2*rgb2; mse2=3*mse(sensorXYZ2,validXYZ) measLab2=xyz2lab(sensorXYZ2',whitePoint'); mseLab2=3*mse(measLab2,validLab) % Third fitting method % Use R G B RG GB BR RR GG BB [colorTrans3] = Model3(validRGB,validXYZ); [rgb3] = Model3Rgb(validRGB); sensorXYZ3=colorTrans3*rgb3; mse3 = 3*mse(sensorXYZ3,validXYZ) measLab3=xyz2lab(sensorXYZ3',whitePoint'); mseLab3=3*mse(measLab3,validLab) % So now the multimillion dollar question is,... % How bad is the error for colors not in the calibration set? % I fear to calibrate with few patches, so I'm going to leave one % patch out each time, fit the model, then calculate the error in % the computed XYZ value for the omitted patch. % Sum up squared XYZ errors in firstOrderXYZError. % measRGB is original full set of checkers. % validCheckers are indices into measRGB. numIters=size(validCheckers,2); indices=[1:numIters]; checkerSubset=indices(2:numIters); computedXYZ1=zeros(size(validXYZ)); computedXYZ2=computedXYZ1; computedXYZ3=computedXYZ1; % i = number of the checker that was left out for i = 1:numIters meterXYZ=validXYZ(:,checkerSubset); cameraRGB=validRGB(:,checkerSubset); % Model 1 colorTrans = cameraRGB'\meterXYZ'; colorTrans = colorTrans'; computedXYZ1(:,i)=colorTrans*validRGB(:,i); % Model 2 cameraRGB2=Model2Rgb(cameraRGB); colorTrans2=Model2(cameraRGB,meterXYZ); computedXYZ2(:,i)=colorTrans2*rgb2(:,i); % Model 3 cameraRGB3=Model3Rgb(cameraRGB); colorTrans3=Model3(cameraRGB,meterXYZ); computedXYZ3(:,i)=colorTrans3*rgb3(:,i); if(i