blackbody | Compute blackbody emission spectrum
%
% E = BLACKBODY(LAMBDA, T) is the blackbody radiation power density (W/m^3)
% at the wavelength LAMBDA (m) and temperature T (K).
%
% If LAMBDA is a column vector, then E is a column vector whose
% elements correspond to to those in LAMBDA. For example:
%
% l = [380:10:700]'*1e-9; % visible spectrum
% e = blackbody(l, 6500); % emission of sun
% plot(l, e)
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function e = blackbody(lam, T)
C1 = 5.9634e-17;
C2 = 1.4387e-2;
lam = lam(:);
e = 2 * C1 ./ (lam.^5 .* (exp(C2/T./lam) - 1)); |
ccdresponse | CCD spectral response
%
% R = CCDRESPONSE(LAMBDA) is the spectral response of a typical silicon
% imaging sensor at the wavelength LAMBDA. The response is normalized
% in the range 0 to 1. If LAMBDA is a vector then R is a vector of the
% same length whose elements are the response at the corresponding element
% of LAMBDA.
%
% Reference::
% Data taken from an ancient Fairchild data book for a sensor with no IR filter fitted.
%
% See also RLUMINOS.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function cc = ccdresponse(lam)
tab = [300 0
350 0
400 5
500 30
600 60
700 85
800 100
900 85
1000 50
];
cc = spline(tab(:,1)*1e-9,tab(:,2),lam)/100; |
cie_primaries | Define CIE primary colors
%
% P = CIE_PRIMARIES() is a 3-vector with the wavelengths (m) of the
% CIE 1976 red, green and blue primaries respectively.
function p = cie_primaries
p = [700, 546.1, 435.8]*1e-9; |
cmfrgb | Color matching function
%
% RGB = CMFRGB(LAMBDA) is the CIE color matching function for illumination
% at wavelength LAMBDA. If LAMBDA is a vector then each row of XYZ
% is the color matching function of the corresponding element of LAMBDA.
%
% RGB = CMFRGB(LAMBDA, E) is the CIE color matching function for an illumination
% spectrum E. E and LAMBDA are vectors of the same length and the elements of E
% represent the intensity of light at the corresponding wavelength in LAMBDA.
%
% Note::
% - the color matching function is the tristimulus required to match a
% particular wavelength excitation.
% - data from http://cvrl.ioo.ucl.ac.uk
% - The Stiles & Burch 2-deg CMFs are based on measurements made on
% 10 observers. The data are referred to as pilot data, but probably
% represent the best estimate of the 2 deg CMFs, since, unlike the CIE
% 2 deg functions (which were reconstructed from chromaticity data),
% they were measured directly.
% - From Table I(5.5.3) of Wyszecki & Stiles (1982). (Table 1(5.5.3)
% of Wyszecki & Stiles (1982) gives the Stiles & Burch functions in
% 250 cm-1 steps, while Table I(5.5.3) of Wyszecki & Stiles (1982)
% is gives them in interpolated 1 nm steps.)
% - These CMFs differ slightly from those of Stiles & Burch (1955). As
% noted in footnote a on p. 335 of Table 1(5.5.3) of Wyszecki &
% Stiles (1982), the CMFs have been "corrected in accordance with
% instructions given by Stiles & Burch (1959)" and renormalized to
% primaries at 15500 (645.16), 19000 (526.32), and 22500 (444.44) cm-1
%
% See also CCXYZ.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function rgb = cmfrgb(lambda, spect)
ciedat = [
390e-9, 1.83970e-003, -4.53930e-004, 1.21520e-002
395e-9, 4.61530e-003, -1.04640e-003, 3.11100e-002
400e-9, 9.62640e-003, -2.16890e-003, 6.23710e-002
405e-9, 1.89790e-002, -4.43040e-003, 1.31610e-001
410e-9, 3.08030e-002, -7.20480e-003, 2.27500e-001
415e-9, 4.24590e-002, -1.25790e-002, 3.58970e-001
420e-9, 5.16620e-002, -1.66510e-002, 5.23960e-001
425e-9, 5.28370e-002, -2.12400e-002, 6.85860e-001
430e-9, 4.42870e-002, -1.99360e-002, 7.96040e-001
435e-9, 3.22200e-002, -1.60970e-002, 8.94590e-001
440e-9, 1.47630e-002, -7.34570e-003, 9.63950e-001
445e-9, -2.33920e-003, 1.36900e-003, 9.98140e-001
450e-9, -2.91300e-002, 1.96100e-002, 9.18750e-001
455e-9, -6.06770e-002, 4.34640e-002, 8.24870e-001
460e-9, -9.62240e-002, 7.09540e-002, 7.85540e-001
465e-9, -1.37590e-001, 1.10220e-001, 6.67230e-001
470e-9, -1.74860e-001, 1.50880e-001, 6.10980e-001
475e-9, -2.12600e-001, 1.97940e-001, 4.88290e-001
480e-9, -2.37800e-001, 2.40420e-001, 3.61950e-001
485e-9, -2.56740e-001, 2.79930e-001, 2.66340e-001
490e-9, -2.77270e-001, 3.33530e-001, 1.95930e-001
495e-9, -2.91250e-001, 4.05210e-001, 1.47300e-001
500e-9, -2.95000e-001, 4.90600e-001, 1.07490e-001
505e-9, -2.97060e-001, 5.96730e-001, 7.67140e-002
510e-9, -2.67590e-001, 7.01840e-001, 5.02480e-002
515e-9, -2.17250e-001, 8.08520e-001, 2.87810e-002
520e-9, -1.47680e-001, 9.10760e-001, 1.33090e-002
525e-9, -3.51840e-002, 9.84820e-001, 2.11700e-003
530e-9, 1.06140e-001, 1.03390e+000, -4.15740e-003
535e-9, 2.59810e-001, 1.05380e+000, -8.30320e-003
540e-9, 4.19760e-001, 1.05120e+000, -1.21910e-002
545e-9, 5.92590e-001, 1.04980e+000, -1.40390e-002
550e-9, 7.90040e-001, 1.03680e+000, -1.46810e-002
555e-9, 1.00780e+000, 9.98260e-001, -1.49470e-002
560e-9, 1.22830e+000, 9.37830e-001, -1.46130e-002
565e-9, 1.47270e+000, 8.80390e-001, -1.37820e-002
570e-9, 1.74760e+000, 8.28350e-001, -1.26500e-002
575e-9, 2.02140e+000, 7.46860e-001, -1.13560e-002
580e-9, 2.27240e+000, 6.49300e-001, -9.93170e-003
585e-9, 2.48960e+000, 5.63170e-001, -8.41480e-003
590e-9, 2.67250e+000, 4.76750e-001, -7.02100e-003
595e-9, 2.80930e+000, 3.84840e-001, -5.74370e-003
600e-9, 2.87170e+000, 3.00690e-001, -4.27430e-003
605e-9, 2.85250e+000, 2.28530e-001, -2.91320e-003
610e-9, 2.76010e+000, 1.65750e-001, -2.26930e-003
615e-9, 2.59890e+000, 1.13730e-001, -1.99660e-003
620e-9, 2.37430e+000, 7.46820e-002, -1.50690e-003
625e-9, 2.10540e+000, 4.65040e-002, -9.38220e-004
630e-9, 1.81450e+000, 2.63330e-002, -5.53160e-004
635e-9, 1.52470e+000, 1.27240e-002, -3.16680e-004
640e-9, 1.25430e+000, 4.50330e-003, -1.43190e-004
645e-9, 1.00760e+000, 9.66110e-005, -4.08310e-006
650e-9, 7.86420e-001, -1.96450e-003, 1.10810e-004
655e-9, 5.96590e-001, -2.63270e-003, 1.91750e-004
660e-9, 4.43200e-001, -2.62620e-003, 2.26560e-004
665e-9, 3.24100e-001, -2.30270e-003, 2.15200e-004
670e-9, 2.34550e-001, -1.87000e-003, 1.63610e-004
675e-9, 1.68840e-001, -1.44240e-003, 9.71640e-005
680e-9, 1.20860e-001, -1.07550e-003, 5.10330e-005
685e-9, 8.58110e-002, -7.90040e-004, 3.52710e-005
690e-9, 6.02600e-002, -5.67650e-004, 3.12110e-005
695e-9, 4.14800e-002, -3.92740e-004, 2.45080e-005
700e-9, 2.81140e-002, -2.62310e-004, 1.65210e-005
705e-9, 1.91170e-002, -1.75120e-004, 1.11240e-005
710e-9, 1.33050e-002, -1.21400e-004, 8.69650e-006
715e-9, 9.40920e-003, -8.57600e-005, 7.43510e-006
720e-9, 6.51770e-003, -5.76770e-005, 6.10570e-006
725e-9, 4.53770e-003, -3.90030e-005, 5.02770e-006
730e-9, 3.17420e-003, -2.65110e-005, 4.12510e-006];
rgb = interp1(ciedat(:,1), ciedat(:,2:4), lambda, 'spline', 0);
if nargin == 2,
rgb = spect(:)' * rgb / numrows(ciedat);
end |
cmfxyz | matching function
%
% XYZ = CMFXYZ(LAMBDA) is the CIE XYZ color matching function for
% illumination at wavelength LAMBDA. If LAMBDA is a vector then each row of
% XYZ is the color matching function of the corresponding element of LAMBDA.
%
% XYZ = CMFXYZ(LAMBDA, E) is the CIE XYZ color matching function for an
% illumination spectrum E. E and LAMBDA are vectors of the same length and
% the elements of E represent the intensity of light at the corresponding
% wavelength in LAMBDA.
%
% Note::
% - the color matching function is the XYZ tristimulus required to match a
% particular wavelength excitation.
% - CIE 1931 2-deg XYZ CMFs from cvrl.ioo.ucl.ac.uk
%
% See also CMFRGB, CCXYZ.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function XYZ = cmfxyz(lambda, spect)
if true
ciedat = [
0 0 0 0
% CIE 1931 2-deg XYZ CMFs from cvrl.ioo.ucl.ac.uk
360,0.000129900000,0.000003917000,0.000606100000
365,0.000232100000,0.000006965000,0.001086000000
370,0.000414900000,0.000012390000,0.001946000000
375,0.000741600000,0.000022020000,0.003486000000
380,0.001368000000,0.000039000000,0.006450001000
385,0.002236000000,0.000064000000,0.010549990000
390,0.004243000000,0.000120000000,0.020050010000
395,0.007650000000,0.000217000000,0.036210000000
400,0.014310000000,0.000396000000,0.067850010000
405,0.023190000000,0.000640000000,0.110200000000
410,0.043510000000,0.001210000000,0.207400000000
415,0.077630000000,0.002180000000,0.371300000000
420,0.134380000000,0.004000000000,0.645600000000
425,0.214770000000,0.007300000000,1.039050100000
430,0.283900000000,0.011600000000,1.385600000000
435,0.328500000000,0.016840000000,1.622960000000
440,0.348280000000,0.023000000000,1.747060000000
445,0.348060000000,0.029800000000,1.782600000000
450,0.336200000000,0.038000000000,1.772110000000
455,0.318700000000,0.048000000000,1.744100000000
460,0.290800000000,0.060000000000,1.669200000000
465,0.251100000000,0.073900000000,1.528100000000
470,0.195360000000,0.090980000000,1.287640000000
475,0.142100000000,0.112600000000,1.041900000000
480,0.095640000000,0.139020000000,0.812950100000
485,0.057950010000,0.169300000000,0.616200000000
490,0.032010000000,0.208020000000,0.465180000000
495,0.014700000000,0.258600000000,0.353300000000
500,0.004900000000,0.323000000000,0.272000000000
505,0.002400000000,0.407300000000,0.212300000000
510,0.009300000000,0.503000000000,0.158200000000
515,0.029100000000,0.608200000000,0.111700000000
520,0.063270000000,0.710000000000,0.078249990000
525,0.109600000000,0.793200000000,0.057250010000
530,0.165500000000,0.862000000000,0.042160000000
535,0.225749900000,0.914850100000,0.029840000000
540,0.290400000000,0.954000000000,0.020300000000
545,0.359700000000,0.980300000000,0.013400000000
550,0.433449900000,0.994950100000,0.008749999000
555,0.512050100000,1.000000000000,0.005749999000
560,0.594500000000,0.995000000000,0.003900000000
565,0.678400000000,0.978600000000,0.002749999000
570,0.762100000000,0.952000000000,0.002100000000
575,0.842500000000,0.915400000000,0.001800000000
580,0.916300000000,0.870000000000,0.001650001000
585,0.978600000000,0.816300000000,0.001400000000
590,1.026300000000,0.757000000000,0.001100000000
595,1.056700000000,0.694900000000,0.001000000000
600,1.062200000000,0.631000000000,0.000800000000
605,1.045600000000,0.566800000000,0.000600000000
610,1.002600000000,0.503000000000,0.000340000000
615,0.938400000000,0.441200000000,0.000240000000
620,0.854449900000,0.381000000000,0.000190000000
625,0.751400000000,0.321000000000,0.000100000000
630,0.642400000000,0.265000000000,0.000049999990
635,0.541900000000,0.217000000000,0.000030000000
640,0.447900000000,0.175000000000,0.000020000000
645,0.360800000000,0.138200000000,0.000010000000
650,0.283500000000,0.107000000000,0.000000000000
655,0.218700000000,0.081600000000,0.000000000000
660,0.164900000000,0.061000000000,0.000000000000
665,0.121200000000,0.044580000000,0.000000000000
670,0.087400000000,0.032000000000,0.000000000000
675,0.063600000000,0.023200000000,0.000000000000
680,0.046770000000,0.017000000000,0.000000000000
685,0.032900000000,0.011920000000,0.000000000000
690,0.022700000000,0.008210000000,0.000000000000
695,0.015840000000,0.005723000000,0.000000000000
700,0.011359160000,0.004102000000,0.000000000000
705,0.008110916000,0.002929000000,0.000000000000
710,0.005790346000,0.002091000000,0.000000000000
715,0.004109457000,0.001484000000,0.000000000000
720,0.002899327000,0.001047000000,0.000000000000
725,0.002049190000,0.000740000000,0.000000000000
730,0.001439971000,0.000520000000,0.000000000000
735,0.000999949300,0.000361100000,0.000000000000
740,0.000690078600,0.000249200000,0.000000000000
745,0.000476021300,0.000171900000,0.000000000000
750,0.000332301100,0.000120000000,0.000000000000
755,0.000234826100,0.000084800000,0.000000000000
760,0.000166150500,0.000060000000,0.000000000000
765,0.000117413000,0.000042400000,0.000000000000
770,0.000083075270,0.000030000000,0.000000000000
775,0.000058706520,0.000021200000,0.000000000000
780,0.000041509940,0.000014990000,0.000000000000
785,0.000029353260,0.000010600000,0.000000000000
790,0.000020673830,0.000007465700,0.000000000000
795,0.000014559770,0.000005257800,0.000000000000
800,0.000010253980,0.000003702900,0.000000000000
805,0.000007221456,0.000002607800,0.000000000000
810,0.000005085868,0.000001836600,0.000000000000
815,0.000003581652,0.000001293400,0.000000000000
820,0.000002522525,0.000000910930,0.000000000000
825,0.000001776509,0.000000641530,0.000000000000
830,0.000001251141,0.000000451810,0.000000000000
1e9 0 0 0];
else
ciedat = [
0 0 0 0
360.000000000000,0.000000122200,0.000000013398,0.000000535027
365.000000000000,0.000000919270,0.000000100650,0.000004028300
370.000000000000,0.000005958600,0.000000651100,0.000026143700
375.000000000000,0.000033266000,0.000003625000,0.000146220000
380.000000000000,0.000159952000,0.000017364000,0.000704776000
385.000000000000,0.000662440000,0.000071560000,0.002927800000
390.000000000000,0.002361600000,0.000253400000,0.010482200000
395.000000000000,0.007242300000,0.000768500000,0.032344000000
400.000000000000,0.019109700000,0.002004400000,0.086010900000
405.000000000000,0.043400000000,0.004509000000,0.197120000000
410.000000000000,0.084736000000,0.008756000000,0.389366000000
415.000000000000,0.140638000000,0.014456000000,0.656760000000
420.000000000000,0.204492000000,0.021391000000,0.972542000000
425.000000000000,0.264737000000,0.029497000000,1.282500000000
430.000000000000,0.314679000000,0.038676000000,1.553480000000
435.000000000000,0.357719000000,0.049602000000,1.798500000000
440.000000000000,0.383734000000,0.062077000000,1.967280000000
445.000000000000,0.386726000000,0.074704000000,2.027300000000
450.000000000000,0.370702000000,0.089456000000,1.994800000000
455.000000000000,0.342957000000,0.106256000000,1.900700000000
460.000000000000,0.302273000000,0.128201000000,1.745370000000
465.000000000000,0.254085000000,0.152761000000,1.554900000000
470.000000000000,0.195618000000,0.185190000000,1.317560000000
475.000000000000,0.132349000000,0.219940000000,1.030200000000
480.000000000000,0.080507000000,0.253589000000,0.772125000000
485.000000000000,0.041072000000,0.297665000000,0.570060000000
490.000000000000,0.016172000000,0.339133000000,0.415254000000
495.000000000000,0.005132000000,0.395379000000,0.302356000000
500.000000000000,0.003816000000,0.460777000000,0.218502000000
505.000000000000,0.015444000000,0.531360000000,0.159249000000
510.000000000000,0.037465000000,0.606741000000,0.112044000000
515.000000000000,0.071358000000,0.685660000000,0.082248000000
520.000000000000,0.117749000000,0.761757000000,0.060709000000
525.000000000000,0.172953000000,0.823330000000,0.043050000000
530.000000000000,0.236491000000,0.875211000000,0.030451000000
535.000000000000,0.304213000000,0.923810000000,0.020584000000
540.000000000000,0.376772000000,0.961988000000,0.013676000000
545.000000000000,0.451584000000,0.982200000000,0.007918000000
550.000000000000,0.529826000000,0.991761000000,0.003988000000
555.000000000000,0.616053000000,0.999110000000,0.001091000000
560.000000000000,0.705224000000,0.997340000000,0.000000000000
565.000000000000,0.793832000000,0.982380000000,0.000000000000
570.000000000000,0.878655000000,0.955552000000,0.000000000000
575.000000000000,0.951162000000,0.915175000000,0.000000000000
580.000000000000,1.014160000000,0.868934000000,0.000000000000
585.000000000000,1.074300000000,0.825623000000,0.000000000000
590.000000000000,1.118520000000,0.777405000000,0.000000000000
595.000000000000,1.134300000000,0.720353000000,0.000000000000
600.000000000000,1.123990000000,0.658341000000,0.000000000000
605.000000000000,1.089100000000,0.593878000000,0.000000000000
610.000000000000,1.030480000000,0.527963000000,0.000000000000
615.000000000000,0.950740000000,0.461834000000,0.000000000000
620.000000000000,0.856297000000,0.398057000000,0.000000000000
625.000000000000,0.754930000000,0.339554000000,0.000000000000
630.000000000000,0.647467000000,0.283493000000,0.000000000000
635.000000000000,0.535110000000,0.228254000000,0.000000000000
640.000000000000,0.431567000000,0.179828000000,0.000000000000
645.000000000000,0.343690000000,0.140211000000,0.000000000000
650.000000000000,0.268329000000,0.107633000000,0.000000000000
655.000000000000,0.204300000000,0.081187000000,0.000000000000
660.000000000000,0.152568000000,0.060281000000,0.000000000000
665.000000000000,0.112210000000,0.044096000000,0.000000000000
670.000000000000,0.081260600000,0.031800400000,0.000000000000
675.000000000000,0.057930000000,0.022601700000,0.000000000000
680.000000000000,0.040850800000,0.015905100000,0.000000000000
685.000000000000,0.028623000000,0.011130300000,0.000000000000
690.000000000000,0.019941300000,0.007748800000,0.000000000000
695.000000000000,0.013842000000,0.005375100000,0.000000000000
700.000000000000,0.009576880000,0.003717740000,0.000000000000
705.000000000000,0.006605200000,0.002564560000,0.000000000000
710.000000000000,0.004552630000,0.001768470000,0.000000000000
715.000000000000,0.003144700000,0.001222390000,0.000000000000
720.000000000000,0.002174960000,0.000846190000,0.000000000000
725.000000000000,0.001505700000,0.000586440000,0.000000000000
730.000000000000,0.001044760000,0.000407410000,0.000000000000
735.000000000000,0.000727450000,0.000284041000,0.000000000000
740.000000000000,0.000508258000,0.000198730000,0.000000000000
745.000000000000,0.000356380000,0.000139550000,0.000000000000
750.000000000000,0.000250969000,0.000098428000,0.000000000000
755.000000000000,0.000177730000,0.000069819000,0.000000000000
760.000000000000,0.000126390000,0.000049737000,0.000000000000
765.000000000000,0.000090151000,0.000035540500,0.000000000000
770.000000000000,0.000064525800,0.000025486000,0.000000000000
775.000000000000,0.000046339000,0.000018338400,0.000000000000
780.000000000000,0.000033411700,0.000013249000,0.000000000000
785.000000000000,0.000024209000,0.000009619600,0.000000000000
790.000000000000,0.000017611500,0.000007012800,0.000000000000
795.000000000000,0.000012855000,0.000005129800,0.000000000000
800.000000000000,0.000009413630,0.000003764730,0.000000000000
805.000000000000,0.000006913000,0.000002770810,0.000000000000
810.000000000000,0.000005093470,0.000002046130,0.000000000000
815.000000000000,0.000003767100,0.000001516770,0.000000000000
820.000000000000,0.000002795310,0.000001128090,0.000000000000
825.000000000000,0.000002082000,0.000000842160,0.000000000000
830.000000000000,0.000001553140,0.000000629700,0.000000000000
1e9 0 0 0];
end
XYZ = interp1(ciedat(:,1)*1e-9, ciedat(:,2:4), lambda, 'pchip', 0);
if nargin == 2,
XYZ = spect(:)' * XYZ;
end |
colordistance | Return distance in RG colorspace
%
% D = COLORDISTANCE(RGB, RG) is the distance on the rg-chromaticity plane from
% coordinate RG=[r,g] to every pixel in the color image RGB. D is an image
% with the same dimensions as RGB and the value of each pixel is the the color
% space distance of the corresponding pixel in RGB.
%
% Note::
% - The output image could be thresholded to determine color similarity.
%
% See also COLORSEG.
% Peter Corke 2005
%
function s = colordistance(im, rg)
% convert image to (r,g) coordinates
rgb = sum(im, 3);
r = im(:,:,1) ./ rgb;
g = im(:,:,2) ./ rgb;
% compute the Euclidean color space distance
d = (r - rg(1)).^2 + (g - rg(2)).^2;
if nargout == 0,
idisp(d)
else
s = d;
end |
colorname | Map between color names and RGB values
%
% RGB = COLORNAME(NAME) is the rgb-tristimulus value corresponding to the
% color specified by the string NAME.
%
% NAME = COLORNAME(RGB) is a string giving the name of the color that is
% closest (Euclidean) to the given rgb-tristimulus value.
%
% Notes::
% - Based on the standard X11 color database rgb.txt.
% - tristimulus values are in the range 0 to 1
function r = colorname(a, opt)
if nargin < 2
opt = '';
end
persistent rgbtable;
% ensure that the database is loaded
if isempty(rgbtable),
% load mapping table from file
fprintf('loading rgb.txt\n');
f = fopen('private/rgb.txt', 'r');
k = 0;
rgb = [];
names = {};
xy = [];
while ~feof(f),
line = fgets(f);
if line(1) == '#',
continue;
end
[A,count,errm,next] = sscanf(line, '%d %d %d');
if count == 3,
k = k + 1;
rgb(k,:) = A' / 255.0;
names{k} = lower( strtrim(line(next:end)) );
xy = xyz2xy( colorspace('RGB->XYZ', rgb) );
end
end
s.rgb = rgb;
s.names = names;
s.xy = xy;
rgbtable = s;
end
if isstr(a)
% map name to rgb
if a(1) == '?'
% just do a wildcard lookup
r = namelookup(rgbtable, a(2:end));
else
r = name2rgb(rgbtable, a, opt);
end
elseif iscell(a)
% map multiple names to rgb
r = [];
for name=a,
rgb = name2rgb(rgbtable, name{1}, opt);
if isempty(rgb)
warning('Color %s not found', name{1});
end
r = [r; rgb];
end
else
if numel(a) == 3
r = rgb2name(rgbtable, a(:)');
elseif numcols(a) == 2 && strcmp(opt, 'xy')
% convert xy to a name
r = {};
for k=1:numrows(a),
r{k} = xy2name(rgbtable, a(k,:));
end
elseif numcols(a) == 3 && ~strcmp(opt, 'xy')
% convert RGB data to a name
r = {};
for k=1:numrows(a),
r{k} = rgb2name(rgbtable, a(k,:));
end
end
end
end
function r = namelookup(table, s)
s = lower(s); % all matching done in lower case
r = {};
count = 1;
for k=1:length(table.names),
if ~isempty( findstr(table.names{k}, s) )
r{count} = table.names{k};
count = count + 1;
end
end
end
function r = name2rgb(table, s, opt)
s = lower(s); % all matching done in lower case
for k=1:length(table.names),
if strcmp(s, table.names(k)),
r = table.rgb(k,:);
if strcmp(opt, 'xy')
r
XYZ = colorspace('RGB->XYZ', r);
XYZ
r = tristim2cc(XYZ);
r
end
return;
end
end
r = [];
end
function r = rgb2name(table, v)
d = table.rgb - ones(numrows(table.rgb),1) * v;
n = colnorm(d');
[z,k] = min(n);
r = table.names{k};
end
function r = xy2name(table, v)
d = table.xy - ones(numrows(table.xy),1) * v;
n = colnorm(d');
[z,k] = min(n);
r = table.names{k};
end |
colorspace | Convert a color image between color representations.
%
% B = COLORSPACE(S,A) converts the color representation of image A
% where S is a string specifying the conversion. S tells the
% source and destination color spaces, S = 'dest<-src', or
% alternatively, S = 'src->dest'.
%
% [B1,B2,B3] = COLORSPACE(S,A) as above but specifies separate output channels.
%
% COLORSPACE(S,A1,A2,A3) as above but specifies separate input channels.
%
%
% Supported color spaces are:
%
% 'RGB' R'G'B' Red Green Blue (ITU-R BT.709 gamma-corrected)
% 'YPbPr' Luma (ITU-R BT.601) + Chroma
% 'YCbCr'/'YCC' Luma + Chroma ("digitized" version of Y'PbPr)
% 'YUV' NTSC PAL Y'UV Luma + Chroma
% 'YIQ' NTSC Y'IQ Luma + Chroma
% 'YDbDr' SECAM Y'DbDr Luma + Chroma
% 'JPEGYCbCr' JPEG-Y'CbCr Luma + Chroma
% 'HSV'/'HSB' Hue Saturation Value/Brightness
% 'HSL'/'HLS'/'HSI' Hue Saturation Luminance/Intensity
% 'XYZ' CIE XYZ
% 'Lab' CIE L*a*b* (CIELAB)
% 'Luv' CIE L*u*v* (CIELUV)
% 'Lch' CIE L*ch (CIELCH)
%
% Notes::
% - All conversions assume 2 degree observer and D65 illuminant.
% - Color space names are case insensitive.
% - When R'G'B' is the source or destination, it can be omitted. For
% example 'yuv<-' is short for 'yuv<-rgb'.
% - MATLAB uses two standard data formats for R'G'B': double data with
% intensities in the range 0 to 1, and uint8 data with integer-valued
% intensities from 0 to 255. As MATLAB's native datatype, double data is
% the natural choice, and the R'G'B' format used by colorspace. However,
% for memory and computational performance, some functions also operate
% with uint8 R'G'B'. Given uint8 R'G'B' color data, colorspace will
% first cast it to double R'G'B' before processing.
% - If A is an Mx3 array, like a colormap, B will also have size Mx3.
%
% Author::
% Pascal Getreuer 2005-2006
function varargout = colorspace(Conversion,varargin)
%%% Input parsing %%%
if nargin < 2, error('Not enough input arguments.'); end
[SrcSpace,DestSpace] = parse(Conversion);
if nargin == 2
Image = varargin{1};
elseif nargin >= 3
Image = cat(3,varargin{:});
else
error('Invalid number of input arguments.');
end
FlipDims = (size(Image,3) == 1);
if FlipDims, Image = permute(Image,[1,3,2]); end
if ~isa(Image,'double'), Image = double(Image)/255; end
if size(Image,3) ~= 3, error('Invalid input size.'); end
SrcT = gettransform(SrcSpace);
DestT = gettransform(DestSpace);
if ~ischar(SrcT) & ~ischar(DestT)
% Both source and destination transforms are affine, so they
% can be composed into one affine operation
T = [DestT(:,1:3)*SrcT(:,1:3),DestT(:,1:3)*SrcT(:,4)+DestT(:,4)];
Temp = zeros(size(Image));
Temp(:,:,1) = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3) + T(10);
Temp(:,:,2) = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3) + T(11);
Temp(:,:,3) = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3) + T(12);
Image = Temp;
elseif ~ischar(DestT)
Image = rgb(Image,SrcSpace);
Temp = zeros(size(Image));
Temp(:,:,1) = DestT(1)*Image(:,:,1) + DestT(4)*Image(:,:,2) + DestT(7)*Image(:,:,3) + DestT(10);
Temp(:,:,2) = DestT(2)*Image(:,:,1) + DestT(5)*Image(:,:,2) + DestT(8)*Image(:,:,3) + DestT(11);
Temp(:,:,3) = DestT(3)*Image(:,:,1) + DestT(6)*Image(:,:,2) + DestT(9)*Image(:,:,3) + DestT(12);
Image = Temp;
else
Image = feval(DestT,Image,SrcSpace);
end
%%% Output format %%%
if nargout > 1
varargout = {Image(:,:,1),Image(:,:,2),Image(:,:,3)};
else
if FlipDims, Image = permute(Image,[1,3,2]); end
varargout = {Image};
end
return;
function [SrcSpace,DestSpace] = parse(Str)
% Parse conversion argument
if isstr(Str)
Str = lower(strrep(strrep(Str,'-',''),' ',''));
k = find(Str == '>');
if length(k) == 1 % Interpret the form 'src->dest'
SrcSpace = Str(1:k-1);
DestSpace = Str(k+1:end);
else
k = find(Str == '<');
if length(k) == 1 % Interpret the form 'dest<-src'
DestSpace = Str(1:k-1);
SrcSpace = Str(k+1:end);
else
error(['Invalid conversion, ''',Str,'''.']);
end
end
SrcSpace = alias(SrcSpace);
DestSpace = alias(DestSpace);
else
SrcSpace = 1; % No source pre-transform
DestSpace = Conversion;
if any(size(Conversion) ~= 3), error('Transformation matrix must be 3x3.'); end
end
return;
function Space = alias(Space)
Space = strrep(Space,'cie','');
if isempty(Space)
Space = 'rgb';
end
switch Space
case {'ycbcr','ycc'}
Space = 'ycbcr';
case {'hsv','hsb'}
Space = 'hsv';
case {'hsl','hsi','hls'}
Space = 'hsl';
case {'rgb','yuv','yiq','ydbdr','ycbcr','jpegycbcr','xyz','lab','luv','lch'}
return;
end
return;
function T = gettransform(Space)
% Get a colorspace transform: either a matrix describing an affine transform,
% or a string referring to a conversion subroutine
switch Space
case 'ypbpr'
T = [0.299,0.587,0.114,0;-0.1687367,-0.331264,0.5,0;0.5,-0.418688,-0.081312,0];
case 'yuv'
% R'G'B' to NTSC/PAL YUV
% Wikipedia: http://en.wikipedia.org/wiki/YUV
T = [0.299,0.587,0.114,0;-0.147,-0.289,0.436,0;0.615,-0.515,-0.100,0];
case 'ydbdr'
% R'G'B' to SECAM YDbDr
% Wikipedia: http://en.wikipedia.org/wiki/YDbDr
T = [0.299,0.587,0.114,0;-0.450,-0.883,1.333,0;-1.333,1.116,0.217,0];
case 'yiq'
% R'G'B' in [0,1] to NTSC YIQ in [0,1];[-0.595716,0.595716];[-0.522591,0.522591];
% Wikipedia: http://en.wikipedia.org/wiki/YIQ
T = [0.299,0.587,0.114,0;0.595716,-0.274453,-0.321263,0;0.211456,-0.522591,0.311135,0];
case 'ycbcr'
% R'G'B' (range [0,1]) to ITU-R BRT.601 (CCIR 601) Y'CbCr
% Wikipedia: http://en.wikipedia.org/wiki/YCbCr
% Poynton, Equation 3, scaling of R'G'B to Y'PbPr conversion
T = [65.481,128.553,24.966,16;-37.797,-74.203,112.0,128;112.0,-93.786,-18.214,128];
case 'jpegycbcr'
% Wikipedia: http://en.wikipedia.org/wiki/YCbCr
T = [0.299,0.587,0.114,0;-0.168736,-0.331264,0.5,0.5;0.5,-0.418688,-0.081312,0.5]*255;
case {'rgb','xyz','hsv','hsl','lab','luv','lch'}
T = Space;
otherwise
error(['Unknown color space, ''',Space,'''.']);
end
return;
function Image = rgb(Image,SrcSpace)
% Convert to Rec. 709 R'G'B' from 'SrcSpace'
switch SrcSpace
case 'rgb'
return;
case 'hsv'
% Convert HSV to R'G'B'
Image = huetorgb((1 - Image(:,:,2)).*Image(:,:,3),Image(:,:,3),Image(:,:,1));
case 'hsl'
% Convert HSL to R'G'B'
L = Image(:,:,3);
Delta = Image(:,:,2).*min(L,1-L);
Image = huetorgb(L-Delta,L+Delta,Image(:,:,1));
case {'xyz','lab','luv','lch'}
% Convert to CIE XYZ
Image = xyz(Image,SrcSpace);
% Convert XYZ to RGB
T = [3.240479,-1.53715,-0.498535;-0.969256,1.875992,0.041556;0.055648,-0.204043,1.057311];
R = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3); % R
G = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3); % G
B = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3); % B
% Desaturate and rescale to constrain resulting RGB values to [0,1]
AddWhite = -min(min(min(R,G),B),0);
Scale = max(max(max(R,G),B)+AddWhite,1);
R = (R + AddWhite)./Scale;
G = (G + AddWhite)./Scale;
B = (B + AddWhite)./Scale;
% Apply gamma correction to convert RGB to Rec. 709 R'G'B'
Image(:,:,1) = gammacorrection(R); % R'
Image(:,:,2) = gammacorrection(G); % G'
Image(:,:,3) = gammacorrection(B); % B'
otherwise % Conversion is through an affine transform
T = gettransform(SrcSpace);
temp = inv(T(:,1:3));
T = [temp,-temp*T(:,4)];
R = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3) + T(10);
G = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3) + T(11);
B = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3) + T(12);
AddWhite = -min(min(min(R,G),B),0);
Scale = max(max(max(R,G),B)+AddWhite,1);
R = (R + AddWhite)./Scale;
G = (G + AddWhite)./Scale;
B = (B + AddWhite)./Scale;
Image(:,:,1) = R;
Image(:,:,2) = G;
Image(:,:,3) = B;
end
% Clip to [0,1]
Image = min(max(Image,0),1);
return;
function Image = xyz(Image,SrcSpace)
% Convert to CIE XYZ from 'SrcSpace'
WhitePoint = [0.950456,1,1.088754];
switch SrcSpace
case 'xyz'
return;
case 'luv'
% Convert CIE L*uv to XYZ
WhitePointU = (4*WhitePoint(1))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
WhitePointV = (9*WhitePoint(2))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
L = Image(:,:,1);
Y = (L + 16)/116;
Y = invf(Y)*WhitePoint(2);
U = Image(:,:,2)./(13*L + 1e-6*(L==0)) + WhitePointU;
V = Image(:,:,3)./(13*L + 1e-6*(L==0)) + WhitePointV;
Image(:,:,1) = -(9*Y.*U)./((U-4).*V - U.*V); % X
Image(:,:,2) = Y; % Y
Image(:,:,3) = (9*Y - (15*V.*Y) - (V.*Image(:,:,1)))./(3*V); % Z
case {'lab','lch'}
Image = lab(Image,SrcSpace);
% Convert CIE L*ab to XYZ
fY = (Image(:,:,1) + 16)/116;
fX = fY + Image(:,:,2)/500;
fZ = fY - Image(:,:,3)/200;
Image(:,:,1) = WhitePoint(1)*invf(fX); % X
Image(:,:,2) = WhitePoint(2)*invf(fY); % Y
Image(:,:,3) = WhitePoint(3)*invf(fZ); % Z
otherwise % Convert from some gamma-corrected space
% Convert to Rec. 701 R'G'B'
Image = rgb(Image,SrcSpace);
% Undo gamma correction
R = invgammacorrection(Image(:,:,1));
G = invgammacorrection(Image(:,:,2));
B = invgammacorrection(Image(:,:,3));
% Convert RGB to XYZ
T = inv([3.240479,-1.53715,-0.498535;-0.969256,1.875992,0.041556;0.055648,-0.204043,1.057311]);
Image(:,:,1) = T(1)*R + T(4)*G + T(7)*B; % X
Image(:,:,2) = T(2)*R + T(5)*G + T(8)*B; % Y
Image(:,:,3) = T(3)*R + T(6)*G + T(9)*B; % Z
end
return;
function Image = hsv(Image,SrcSpace)
% Convert to HSV
Image = rgb(Image,SrcSpace);
V = max(Image,[],3);
S = (V - min(Image,[],3))./(V + (V == 0));
Image(:,:,1) = rgbtohue(Image);
Image(:,:,2) = S;
Image(:,:,3) = V;
return;
function Image = hsl(Image,SrcSpace)
% Convert to HSL
switch SrcSpace
case 'hsv'
% Convert HSV to HSL
MaxVal = Image(:,:,3);
MinVal = (1 - Image(:,:,2)).*MaxVal;
L = 0.5*(MaxVal + MinVal);
temp = min(L,1-L);
Image(:,:,2) = 0.5*(MaxVal - MinVal)./(temp + (temp == 0));
Image(:,:,3) = L;
otherwise
Image = rgb(Image,SrcSpace); % Convert to Rec. 701 R'G'B'
% Convert R'G'B' to HSL
MinVal = min(Image,[],3);
MaxVal = max(Image,[],3);
L = 0.5*(MaxVal + MinVal);
temp = min(L,1-L);
S = 0.5*(MaxVal - MinVal)./(temp + (temp == 0));
Image(:,:,1) = rgbtohue(Image);
Image(:,:,2) = S;
Image(:,:,3) = L;
end
return;
function Image = lab(Image,SrcSpace)
% Convert to CIE L*a*b* (CIELAB)
WhitePoint = [0.950456,1,1.088754];
switch SrcSpace
case 'lab'
return;
case 'lch'
% Convert CIE L*CH to CIE L*ab
C = Image(:,:,2);
Image(:,:,2) = cos(Image(:,:,3)*pi/180).*C; % a*
Image(:,:,3) = sin(Image(:,:,3)*pi/180).*C; % b*
otherwise
Image = xyz(Image,SrcSpace); % Convert to XYZ
% Convert XYZ to CIE L*a*b*
X = Image(:,:,1)/WhitePoint(1);
Y = Image(:,:,2)/WhitePoint(2);
Z = Image(:,:,3)/WhitePoint(3);
fX = f(X);
fY = f(Y);
fZ = f(Z);
Image(:,:,1) = 116*fY - 16; % L*
Image(:,:,2) = 500*(fX - fY); % a*
Image(:,:,3) = 200*(fY - fZ); % b*
end
return;
function Image = luv(Image,SrcSpace)
% Convert to CIE L*u*v* (CIELUV)
WhitePoint = [0.950456,1,1.088754];
WhitePointU = (4*WhitePoint(1))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
WhitePointV = (9*WhitePoint(2))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
Image = xyz(Image,SrcSpace); % Convert to XYZ
U = (4*Image(:,:,1))./(Image(:,:,1) + 15*Image(:,:,2) + 3*Image(:,:,3));
V = (9*Image(:,:,2))./(Image(:,:,1) + 15*Image(:,:,2) + 3*Image(:,:,3));
Y = Image(:,:,2)/WhitePoint(2);
L = 116*f(Y) - 16;
Image(:,:,1) = L; % L*
Image(:,:,2) = 13*L.*(U - WhitePointU); % u*
Image(:,:,3) = 13*L.*(V - WhitePointV); % v*
return;
function Image = lch(Image,SrcSpace)
% Convert to CIE L*ch
Image = lab(Image,SrcSpace); % Convert to CIE L*ab
H = atan2(Image(:,:,3),Image(:,:,2));
H = H*180/pi + 360*(H < 0);
Image(:,:,2) = sqrt(Image(:,:,2).^2 + Image(:,:,3).^2); % C
Image(:,:,3) = H; % H
return;
function Image = huetorgb(m0,m2,H)
% Convert HSV or HSL hue to RGB
N = size(H);
H = min(max(H(:),0),360)/60;
m0 = m0(:);
m2 = m2(:);
F = H - round(H/2)*2;
M = [m0, m0 + (m2-m0).*abs(F), m2];
Num = length(m0);
j = [2 1 0;1 2 0;0 2 1;0 1 2;1 0 2;2 0 1;2 1 0]*Num;
k = floor(H) + 1;
Image = reshape([M(j(k,1)+(1:Num).'),M(j(k,2)+(1:Num).'),M(j(k,3)+(1:Num).')],[N,3]);
return;
function H = rgbtohue(Image)
% Convert RGB to HSV or HSL hue
[M,i] = sort(Image,3);
i = i(:,:,3);
Delta = M(:,:,3) - M(:,:,1);
Delta = Delta + (Delta == 0);
R = Image(:,:,1);
G = Image(:,:,2);
B = Image(:,:,3);
H = zeros(size(R));
k = (i == 1);
H(k) = (G(k) - B(k))./Delta(k);
k = (i == 2);
H(k) = 2 + (B(k) - R(k))./Delta(k);
k = (i == 3);
H(k) = 4 + (R(k) - G(k))./Delta(k);
H = 60*H + 360*(H < 0);
H(Delta == 0) = nan;
return;
function Rp = gammacorrection(R)
Rp = real(1.099*R.^0.45 - 0.099);
i = (R < 0.018);
Rp(i) = 4.5138*R(i);
return;
function R = invgammacorrection(Rp)
R = real(((Rp + 0.099)/1.099).^(1/0.45));
i = (R < 0.018);
R(i) = Rp(i)/4.5138;
return;
function fY = f(Y)
fY = real(Y.^(1/3));
i = (Y < 0.008856);
fY(i) = Y(i)*(841/108) + (4/29);
return;
function Y = invf(fY)
Y = fY.^3;
i = (Y < 0.008856);
Y(i) = (fY(i) - 4/29)*(108/841);
return; |
lambda2rg | RGG chromaticity coordinates
%
% RGB = CCRGB(LAMBDA) is the rg-chromaticity coordinate for illumination
% at the specific wavelength LAMBDA.
%
% RGB = CCRGB(LAMBDA, E) is the rg-chromaticity coordinate for an illumination
% spectrum E. E and LAMBDA are vectors of the same length and the elements of E
% represent the intensity of light at the corresponding wavelength in LAMBDA.
%
% See also CMFRGB, CCXYZ.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function [r,g] = lambda2rg(lambda, e)
if nargin == 1,
RGB = cmfrgb(lambda);
elseif nargin == 2,
RGB = cmfrgb(lambda, e);
end
cc = tristim2cc(RGB);
if nargout == 1
r = cc;
elseif nargout == 2
r = cc(:,1);
g = cc(:,2);
end |
lambda2xy | [x,y] = lambda2xy(lambda, varargin)
cmf = cmfxyz(lambda, varargin{:});
xy = xyz2xy(cmf);
if nargout == 2
x = xy(1);
y = xy(2);
else
x = xy;
end |
loadspectrum | Load spectrum data
%
% S = LOADSPECTRUM(LAMBDA, FILENAME) is spectral data from file FILENAME
% interpolated to wavelengths specified in LAMBDA (Nx1). S is NxD and the
% elements correspond to spectral data at the corresponding element of LAMBDA.
%
% [S,LAMBDA] = LOADSPECTRUM(LAMBDA, FILENAME) as above but also returns the
% passed wavelength LAMBDA.
%
% Notes::
% - The File is assumed to have its first column as wavelength in metres, the
% remainding columns are linearly interpolated and returned as columns of S.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function [s,lam] = loadspectrum(lambda, filename)
lambda = lambda(:);
tab = load(filename);
s = interp1(tab(:,1), tab(:,2:end), lambda);
if nargout == 2
lam = lambda;
end |
luminos | Photopic luminosity function
%
% P = LUMINOS(LAMBDA) is the photopic luminosity function for the wavelengths
% in LAMBDA. If LAMBDA is a vector, then P is a vector whose elements are the
% luminosity at the corresponding elements of LAMBDA.
%
% Luminosity has units of lumens which are related to the intensity with
% which wavelengths as perceived by the light-adapted human eye.
%
% See also RLUMINOS.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function lu = luminos(lam)
tab = [
3.8000000e-007 0.0000000e+000
3.8500000e-007 1.0000000e-004
3.9000000e-007 1.0000000e-004
3.9500000e-007 2.0000000e-004
4.0000000e-007 4.0000000e-004
4.0500000e-007 6.0000000e-004
4.1000000e-007 1.2000000e-003
4.1500000e-007 2.2000000e-003
4.2000000e-007 4.0000000e-003
4.2500000e-007 7.3000000e-003
4.3000000e-007 1.1600000e-002
4.3500000e-007 1.6800000e-002
4.4000000e-007 2.3000000e-002
4.4500000e-007 2.9800000e-002
4.5000000e-007 3.8000000e-002
4.5500000e-007 4.8000000e-002
4.6000000e-007 6.0000000e-002
4.6500000e-007 7.3900000e-002
4.7000000e-007 9.1000000e-002
4.7500000e-007 1.1260000e-001
4.8000000e-007 1.3900000e-001
4.8500000e-007 1.6930000e-001
4.9000000e-007 2.0800000e-001
4.9500000e-007 2.5860000e-001
5.0000000e-007 3.2300000e-001
5.0500000e-007 4.0730000e-001
5.1000000e-007 5.0300000e-001
5.1500000e-007 6.0820000e-001
5.2000000e-007 7.1000000e-001
5.2500000e-007 7.9320000e-001
5.3000000e-007 8.6200000e-001
5.3500000e-007 9.1490000e-001
5.4000000e-007 9.5400000e-001
5.4500000e-007 9.8030000e-001
5.5000000e-007 9.9500000e-001
5.5500000e-007 1.0002000e+000
5.6000000e-007 9.9500000e-001
5.6500000e-007 9.7860000e-001
5.7000000e-007 9.5200000e-001
5.7500000e-007 9.1540000e-001
5.8000000e-007 8.7000000e-001
5.8500000e-007 8.1630000e-001
5.9000000e-007 7.5700000e-001
5.9500000e-007 6.9490000e-001
6.0000000e-007 6.3100000e-001
6.0500000e-007 5.6680000e-001
6.1000000e-007 5.0300000e-001
6.1500000e-007 4.4120000e-001
6.2000000e-007 3.8100000e-001
6.2500000e-007 3.2100000e-001
6.3000000e-007 2.6500000e-001
6.3500000e-007 2.1700000e-001
6.4000000e-007 1.7500000e-001
6.4500000e-007 1.3820000e-001
6.5000000e-007 1.0700000e-001
6.5500000e-007 8.1600000e-002
6.6000000e-007 6.1000000e-002
6.6500000e-007 4.4600000e-002
6.7000000e-007 3.2000000e-002
6.7500000e-007 2.3200000e-002
6.8000000e-007 1.7000000e-002
6.8500000e-007 1.1900000e-002
6.9000000e-007 8.2000000e-003
6.9500000e-007 5.7000000e-003
7.0000000e-007 4.1000000e-003
7.0500000e-007 2.9000000e-003
7.1000000e-007 2.1000000e-003
7.1500000e-007 1.5000000e-003
7.2000000e-007 1.0000000e-003
7.2500000e-007 7.0000000e-004
7.3000000e-007 5.0000000e-004
7.3500000e-007 4.0000000e-004
7.4000000e-007 3.0000000e-004
7.4500000e-007 2.0000000e-004
7.5000000e-007 1.0000000e-004
7.5500000e-007 1.0000000e-004
7.6000000e-007 1.0000000e-004
7.6500000e-007 0.0000000e+000
7.7000000e-007 0.0000000e+000];
lu = [];
for lmd = lam(:)',
if (lmd < 380e-9) | (lmd > 770e-9)
lu = [lu; 0.0];
else
lu = [lu; 683*interp1(tab(:,1), tab(:,2), lmd)];
end
end; |
rg_addticks | Add wavelength ticks to spectral locus
%
% RG_ADDTICKS() adds wavelength ticks to the spectral locus.
%
% See also XYCOLOURSPACE.
function rg_addticks(lam1, lam2, lamd)
% well spaced points around the locus
lambda = [460:10:540];
lambda = [lambda 560:20:600];
rgb = cmfrgb(lambda*1e-9);
r = rgb(:,1)./sum(rgb')';
g = rgb(:,2)./sum(rgb')';
hold on
plot(r,g, 'o')
hold off
for i=1:numcols(lambda)
text(r(i), g(i), sprintf(' %d', lambda(i)));
end |
rluminos | Relative photopic luminosity function
%
% P = RUMINOS(LAMBDA) is the relative photopic luminosity function for the
% wavelengths in LAMBDA. If LAMBDA is a vector, then P is a vector whose
% elements are the relative luminosity at the corresponding elements of LAMBDA.
%
% Relative luminosity lies in the interval 0 to 1 which indicate the intensity
% with which wavelengths as perceived by the light-adapted human eye.
%
% See also LUMINOS.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function lu = rluminos(lambda)
xyz = cmfxyz(lambda);
lu = xyz(:,2); % photopic luminosity is the Y color matching function |
tristim2cc | Tristimulus to chromaticity coordinates
%
% CC = TRISTIM2CC(TRI) is the chromaticity coordinates corresponding to the
% tristimulus TRI. If TRI is RGB then CC is rg, if TRI is XYZ then CC is xy.
% Multiple tristimulus values can be given as rows of TRI (Nx3) in which case
% the chromaticity coordinates are the corresponding rows of CC (Nx2).
%
% [C1,C2] = TRISTIM2CC(TRI) as above but the chromaticity coordinates are
% provided in separate vectors, each Nx1.
%
% OUT = TRISTIM2CC(IM) is the chromaticity coordinates corresponding to every
% pixel in the tristimulus image IM (HxWx3). OUT (HxWx2) has planes
% corresponding to r and g, or x and y.
%
% [O1,O2] = TRISTIM2CC(IM) as above but the chromaticity images are provided
% as images (HxW).
function [a,b] = tristim2cc(tri)
if ndims(tri) < 3
% each row is R G B or X Y Z
s = sum(tri')';
cc = tri(:,1:2) ./ [s s];
if nargout == 1
a = cc;
else
a = cc(:,1);
b = cc(:,2);
end
else
s = sum(tri, 3);
if nargout == 1
a = tri(:,:,1:2) ./ cat(3, s, s);
else
a = tri(:,:,1) ./ s;
b = tri(:,:,2) ./ s;
end
end |
xycolorspace | Display spectral locus
%
% XYCOLORSPACE() display the spectral locus in terms of CIE x and y
% coordinates.
%
% XYCOLORSPACE(P) as above but plot the points whose xy-chromaticity
% is given by the columns of P.
%
% [IM,AX,AY] = XYCOLORSPACE() as above returns the spectral locus as an
% image IM, with corresponding x- and y-axis coordinates AX and AY
% respectively.
%
% See also RG_ADDTICKS.
% Based on code by Pascal Getreuer 2006
% Demo for colorspace.m - the CIE xyY "tongue"
function [im,ax,ay] = xycolorspace(cxy)
N = 160;
Nx = round(N*0.8);
Ny = round(N*0.9);
e = 0.01;
% Generate colors in the xyY color space
x = linspace(e,0.8-e,Nx);
y = linspace(e,0.9-e,Ny);
[xx,yy] = meshgrid(x,y);
iyy = 1./(yy + 1e-5*(yy == 0));
% Convert from xyY to XYZ
Y = ones(Ny,Nx);
X = Y.*xx.*iyy;
Z = Y.*(1-xx-yy).*iyy;
% Convert from XYZ to R'G'B'
CIEColor = colorspace('rgb<-xyz',cat(3,X,Y,Z));
% define the boundary
lambda = [400:20:700]*1e-9';
xyz = ccxyz(lambda);
xy = xyz(:,1:2);
% Make a smooth boundary with spline interpolation
xi = [interp1(xy(:,1),1:0.25:size(xy,1),'spline'),xy(1,1)];
yi = [interp1(xy(:,2),1:0.25:size(xy,1),'spline'),xy(1,2)];
% create a mask image, colors within the boundary
in = inpolygon(xx, yy, xi,yi);
CIEColor(~cat(3,in,in,in)) = 1; % set outside pixels to white
if nargout == 0
% Render the colors on the tongue
image(x,y,CIEColor)
if nargin == 1
plot_point(cxy, 'k*', 'textsize', 10, 'sequence', 'textcolor', 'k');
end
set(gca, 'Ydir', 'normal');
hold on
plot(xi,yi,'k-');
% for lambda = 400:20:700,
% xyz = ccxyz(lambda*1e-9);
% plot(xyz(1), xyz(2), 'Marker', 'o', 'MarkerFaceColor', 'k', 'MarkerEdgeColor', 'k');
% end
hold off
axis([0,0.8,0,0.9]);
xlabel('x');
ylabel('y');
grid
shg;
else
ax = x;
ay = y;
im = CIEColor;
if nargin == 1
markfeatures(cxy, 0, '*', {10, 'k'});
end
end |
CentralCamera | Central perspective camera class
%
% C = camera default camera, 1024x1024, f=8mm, 10um pixels, camera at
% origin, optical axis is z-axis, u||x, v||y).
% C = camera(f, s, pp, npix, name)
% C = camera(0) f=sx=sy=1, u0=v0=0: normalized coordinates
%
% This camera model assumes central projection, that is, the focal point
% is at z=0 and the image plane is at z=f. The image is not inverted.
%
% The camera coordinate system is:
%
% 0------------> u X
% |
% |
% | + (principal point)
% |
% |
% v Y
% Z-axis is into the page.
%
% Object properties (read/write)
%
% C.f intrinsic: focal length
% C.s intrinsic: pixel size 2x1
% C.pp intrinsic: principal point 2x1
% C.np size of the virtual image plane (pixels) 2x1
%
% C.Tcam extrinsic: pose of the camera
% C.name name of the camera, used for graphical display
%
% Object properties (read only)
%
% C.su pixel width
% C.sv pixel height
% C.u0 principal point, u coordinate
% C.v0 principal point, v coordinate
%
% Object methods
%
% C.fov return camera half field-of-view angles (2x1 rads)
% C.K return camera intrinsic matrix (3x3)
% C.P return camera project matrix for camera pose (3x4)
% C.P(T) return camera intrinsic matrix for specified camera pose (3x4)
%
% C.H(T, n, d) return homography for plane: normal n, distance d (3x3)
% C.F(T) return fundamental matrix from pose 3x3
% C.E(T) return essential matrix from pose 3x3
% C.E(F) return essential matrix from fundamental matrix 3x3
%
% C.plot_epiline(F, p)
%
% C.invH(H) return solutions for camera motion and plane normal
% C.invE(E) return solutions for pose from essential matrix 4x4x4
%
% C.rpy(r,p,y) set camera rpy angles
% C.rpy(rpy)
%
% uv = C.project(P) return image coordinates for world points P
% uv = C.project(P, T) return image coordinates for world points P
% transformed by T prior to projection
%
% P is a list of 3D world points and the corresponding image plane points are
% returned in uv. Each point is represented by a column in P and uv.
%
% If P has 3 columns it is treated as a number of 3D points in world
% coordinates, one point per row.
%
% If POINTS has 6 columns, each row is treated as the start and end 3D
% coordinate for a line segment, in world coordinates.
%
% The optional arguments, T, represents a transformation that can be applied
% to the object data, P, prior to 'imaging'. The camera pose, C.Tcam, is also
% taken into account.
%
% uv = C.plot(P) display image coordinates for world points P
% uv = C.plot(P, T) isplay image coordinates for world points P transformed by T
%
% Points are displayed as a round marker. Lines are displayed as line segments.
% Optionally returns image plane coordinates uv.
%
% C.show
% C.show(name)
%
% Create a graphical camera with name, and pixel dimensions given by C.npix.
% Automatically called on first call to plot().
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
% TODO:
% make a parent imaging class and subclass perspective, fisheye, panocam
% test for points in front of camera and set to NaN if not
% test for points off the image plane and set to NaN if not
% make clipping/test flags
classdef CentralCamera < Camera
properties
f % focal length
k % radial distortion vector
p % tangential distortion parameters
distortion % camera distortion [k1 k2 k3 p1 p2]
end
properties (SetAccess = private)
end
properties (GetAccess = private, SetAccess = private)
end
properties (GetAccess = private)
end
properties (Dependent = true, SetAccess = private)
end
methods
%
% Return a camera intrinsic parameter structure:
% focal length 8mm
% pixel size 10um square
% image size 1024 x 1024
% principal point (512, 512)
function c = CentralCamera(varargin)
% invoke the superclass constructor
c = c@Camera(varargin{:});
c.type = 'central-perspective';
c.perspective = true;
if nargin == 0,
c.name = 'canonic';
% default values
c.f = 1;
c.distortion = [];
elseif nargin == 1 && isa(varargin{1}, 'CentralCamera')
% copy constructor
old = varargin{1};
for p=properties(c)'
% copy each property across, exceptions occur
% for those with protected SetAccess
p = p{1};
try
c = setfield(c, p, getfield(old, p));
end
end
end
if isempty(c.pp) && ~isempty(c.npix)
c.pp = c.npix/2;
elseif isempty(c.pp)
c.pp =[0 0];
end
end
function n = paramSet(c, args)
n = 0;
switch lower(args{1})
case 'focal'
c.f = args{2}; n = 1;
case 'distortion'
v = args{2}; n = 1;
if length(v) ~= 5
error('distortion vector is [k1 k2 k3 p1 p2]');
end
c.distortion = v;
case 'distortion-bouget'
v = args{2}; n = 1;
if length(v) ~= 5
error('distortion vector is [k1 k2 p1 p2 k3]');
end
c.distortion = [v(1) v(2) v(5) v(3) v(4)];;
case 'default'
c.f = 8e-3; % f
c.rho = [10e-6, 10e-6]; % square pixels 10um side
c.npix = [1024, 1024]; % 1Mpix image plane
c.pp = [512, 512]; % principal point in the middle
c.limits = [0 1024 0 1024];
otherwise
error('VisionToolbox:CentralCamera:UnknownOption', 'Unknown option %s', args{1});
end
end
function s = char(c)
s = sprintf('name: %s [%s]', c.name, c.type);
s = strvcat(s, sprintf( ' focal length: %-11.4g', c.f));
if ~isempty(c.distortion)
s = strvcat(s, sprintf(' distortion: k=(%.4g, %.4g, %.4g), p=(%.4g, %.4g)', c.distortion));
end
s = strvcat(s, char@Camera(c) );
end
% intrinsic parameter matrix
function v = K(c)
v = [ c.f/c.rho(1) 0 c.pp(1)
0 c.f/c.rho(2) c.pp(2)
0 0 1%/c.f
] ;
end
% camera calibration or projection matrix
function v = C(c, Tcam)
if nargin == 1,
Tcam = c.T;
end
if isempty(c.rho)
rho = [1 1];
else
rho = c.rho;
end
v = [ c.f/rho(1) 0 c.pp(1) 0
0 c.f/rho(2) c.pp(2) 0
0 0 1 0
] * inv(Tcam);
end
function HH = H(c, T, n, d)
if (d < 0) || (n(3) < 0)
error('plane distance must be > 0 and normal away from camera');
end
% T is the transform from view 1 to view 2
% (R,t) is the inverse
T = inv(T);
R = t2r( T );
t = transl( T );
HH = R + 1.0/d*t*n(:)';
% now apply the camera intrinsics
K = c.K
HH = K * HH * inv(K);
HH = HH / HH(3,3); % normalize it
end
function s = invH(c, H, varargin)
if nargout == 0
invhomog(H, 'K', c.K, varargin{:});
else
s = invhomog(H, 'K', c.K, varargin{:});
end
end
function fmatrix = F(c, X)
% cam.F(X)
% cam.F(cam2) with
% REF: An Invitation to 3D geometry, p.177
% T is the pose for view 1
% c.T is the pose for view 2
if ishomog(X)
E = c.E(X);
K = c.K();
fmatrix = inv(K)' * E * inv(K);
elseif isa(X, 'Camera')
% use relative pose and camera parameters of
E = c.E(X);
K1 = c.K;
K2 = X.K();
fmatrix = inv(K2)' * E * inv(K1);
end
end
% E = cam1.E(F) convert F matrix to E using camera intrinsics
% E = cam1.E(T12) compute E for second view displaced by T12 from first
% E = cam1.E(cam2) compute E for second view given by cam2 and first view given
% by cam1
function ematrix = E(c, X)
% essential matrix from pose. Assume the first view is associated
% with the passed argument, either a hom.trans or a camera.
% The second view is Tcam of this object.
if ismatrix(X) && all(size(X) == [3 3]),
% essential matrix from F matrix
F = X;
K = c.K();
ematrix = K'*F*K;
return;
elseif isa(X, 'Camera')
T21 = inv(X.T) * c.T;
elseif ishomog(X)
T21 = inv(X);
else
error('unknown argument type');
end
% T = (R,t) is the transform from transform from 2 --> 1
%
% as per Hartley & Zisserman p.244
% P' = K' [R t] where [R t] is inverse camera pose, PIC book (8.5)
%
% as per Ma etal. p.100
% X2 = RX1 + t, so [R t] is the transform from C2 to C1
[R,t] = tr2rt(T21);
% REF: An Invitation to 3D geometry, p.177
% E = Sk(t) R
ematrix = skew(t) * R;
end
function s = invE(c, E, P)
% REF: Hartley & Zisserman, Chap 9, p. 259
% we return T from view 1 to view 2
[U,S,V] = svd(E);
% singular values are (sigma, sigma, 0)
if 0
% H&Z solution
W = [0 -1 0; 1 0 0; 0 0 1]; % rotz(pi/2)
t = U(:,3);
R1 = U*W*V';
if det(R1) < 0,
disp('flip');
V = -V;
R1 = U*W*V';
det(R1)
end
R2 = U*W'*V';
% we need to invert the solutions since our definition of pose is
% from initial camera to the final camera
s(:,:,1) = inv([R1 t; 0 0 0 1]);
s(:,:,2) = inv([R1 -t; 0 0 0 1]);
s(:,:,3) = inv([R2 t; 0 0 0 1]);
s(:,:,4) = inv([R2 -t; 0 0 0 1]);
else
% Ma etal solution, p116, p120-122
% Fig 5.2 (p113), is wrong, (R,t) is from camera 2 to 1
if det(V) < 0
V = -V;
S = -S;
end
if det(U) < 0
U = -U;
S = -S;
end
R1 = U*rotz(pi/2)'*V';
R2 = U*rotz(-pi/2)'*V';
t1 = vex(U*rotz(pi/2)*S*U');
t2 = vex(U*rotz(-pi/2)*S*U');
% invert (R,t) so its from camera 1 to 2
s(:,:,1) = inv( [R1 t1; 0 0 0 1] );
s(:,:,2) = inv( [R2 t2; 0 0 0 1] );
end
if nargin > 2
for i=1:size(s,3)
if ~any(isnan(c.project(P, 'Tcam', s(:,:,i))))
s = s(:,:,i);
fprintf('solution %d is good\n', i);
return;
end
end
warning('no solution has given point in front of camera');
end
end
% plot a line specified in theta-rho format
function plot_line_tr(cam, lines, varargin)
x = get(cam.h_image, 'XLim');
y = get(cam.h_image, 'YLim');
% plot it
for i=1:numcols(lines)
theta = lines(1,i);
rho = lines(2,i);
%fprintf('theta = %f, d = %f\n', line.theta, line.rho);
if abs(cos(theta)) > 0.5,
% horizontalish lines
plot(x, -x*tan(theta) + rho/cos(theta), varargin{:}, 'Parent', cam.h_image);
else
% verticalish lines
plot( -y/tan(theta) + rho/sin(theta), y, varargin{:}, 'Parent', cam.h_image);
end
end
end
function handles = plot_epiline(c, F, p, varargin)
% for all input points
l = F * e2h(p);
c.line(l, varargin{:});
end
% return field-of-view angle for x and y direction (rad)
function th = fov(c)
th = 2*atan(c.npix/2.*c.rho / c.f);
end
function uv = project(c, P, varargin)
% do the camera perspective transform
% P is 3xN matrix of points to plot
opt.Tobj = [];
opt.Tcam = [];
[opt,arglist] = tb_optparse(opt, varargin);
np = numcols(P);
if isempty(opt.Tcam)
opt.Tcam = c.T;
end
if ndims(opt.Tobj) == 3 && ndims(opt.Tcam) == 3
error('cannot animate object and camera simultaneously');
end
if ndims(opt.Tobj) == 3
% animate object motion, static camera
% get camera matrix for this camera pose
C = c.C(opt.Tcam);
% make the world points homogeneous
if numrows(P) == 3
P = e2h(P);
end
for frame=1:size(opt.Tobj,3)
% transform all the points to camera frame
X = C * opt.Tobj(:,:,frame) * P; % project them
X(3,X(3,:)<0) = NaN; % points behind the camera are set to NaN
X = h2e(X); % convert to Euclidean coordinates
if c.noise
% add Gaussian noise with specified standard deviation
X = X + diag(c.noise) * randn(size(X));
end
uv(:,:,frame) = X;
end
else
% animate camera, static object
% transform the object
if ~isempty(opt.Tobj)
P = transformp(opt.Tobj, P);
end
% make the world points homogeneous
if numrows(P) == 3
P = e2h(P);
end
for frame=1:size(opt.Tcam,3)
C = c.C(opt.Tcam(:,:,frame));
% transform all the points to camera frame
X = C * P; % project them
X(3,X(3,:)<0) = NaN; % points behind the camera are set to NaN
X = h2e(X); % convert to Euclidean coordinates
if c.noise
% add Gaussian noise with specified standard deviation
X = X + diag(c.noise) * randn(size(X));
end
uv(:,:,frame) = X;
end
end
end
function r = ray(cam, p)
% from HZ p 162
C = cam.C();
Mi = inv(C(1:3,1:3));
p4 = C(:,4);
for i=1:numcols(p)
r(i) = Ray3D(-Mi*p4, Mi*e2h(p(:,i)));
end
end
function hg = drawCamera(cam, s, varargin)
hold on
if nargin == 0
s = 1;
end
s = s/3;
opt.color = 'b';
opt.mode = {'solid', 'mesh'};
opt.label = false;
opt = tb_optparse(opt, varargin);
% create a new transform group
hg = hgtransform;
% the box is centred at the origin and its centerline parallel to the
% z-axis. Its z-extent is -bh/2 to bh/2.
bw = 0.5; % half width of the box
bh = 1.2; % height of the box
cr = 0.4; % cylinder radius
ch = 0.8; % cylinder height
cn = 16; % number of facets of cylinder
a = 3; % length of axis line segments
opt.scale = s;
opt.parent = hg;
% draw the box part of the camera
r = bw*[1; 1];
x = r * [1 1 -1 -1 1];
y = r * [1 -1 -1 1 1];
z = [-bh; bh]/2 * ones(1,5);
draw(x,y,z, opt);
% draw top and bottom of box
x = bw * [-1 1; -1 1];
y = bw * [1 1; -1 -1];
z = [1 1; 1 1];
draw(x,y,-bh/2*z, opt);
draw(x,y,bh/2*z, opt);
% draw the lens
[x,y,z] = cylinder(cr, cn);
z = ch*z+bh/2;
h = draw(x,y,z, opt);
set(h, 'BackFaceLighting', 'unlit');
% draw the x-, y- and z-axes
h = plot3([0,a*s], [0,0], [0,0], 'k')
set(h, 'Parent', hg);
h = plot3([0,0], [0,a*s], [0,0], 'k')
set(h, 'Parent', hg);
h = plot3([0,0], [0,0], [0,a*s], 'k')
set(h, 'Parent', hg);
if opt.label
h = text( a*s,0,0, cam.name);
set(h, 'Parent', hg);
end
hold off
function h = draw(x, y, z, opt)
s = opt.scale;
switch opt.mode
case 'solid'
h = surf(x*s,y*s,z*s, 'FaceColor', opt.color);
case 'surfl'
h = surfl(x*s,y*s,z*s, 'FaceColor', opt.color);
case 'mesh'
h = mesh(x*s,y*s,z*s, 'EdgeColor', opt.color);
end
set(h, 'Parent', opt.parent);
end
end
end % methods
end % class |
CatadioptricCamera | Catadioptric camera class
%
% C = camera default camera, 1024x1024, f=8mm, 10um pixels, camera at
% origin, optical axis is z-axis, u||x, v||y).
% C = camera(f, s, pp, npix, name)
% C = camera(0) f=sx=sy=1, u0=v0=0: normalized coordinates
%
% This camera model assumes central projection, that is, the focal point
% is at z=0 and the image plane is at z=f. The image is not inverted.
%
% The camera coordinate system is:
%
% 0------------> u X
% |
% |
% | + (principal point)
% |
% |
% v Y
% Z-axis is into the page.
%
% Object properties (read/write)
%
% C.f intrinsic: focal length
% C.s intrinsic: pixel size 2x1
% C.pp intrinsic: principal point 2x1
% C.np size of the virtual image plane (pixels) 2x1
%
% C.Tcam extrinsic: pose of the camera
% C.name name of the camera, used for graphical display
%
% Object properties (read only)
%
% C.su pixel width
% C.sv pixel height
% C.u0 principal point, u coordinate
% C.v0 principal point, v coordinate
%
% Object methods
%
% C.fov return camera half field-of-view angles (2x1 rads)
% C.K return camera intrinsic matrix (3x3)
% C.P return camera project matrix for camera pose (3x4)
% C.P(T) return camera intrinsic matrix for specified camera pose (3x4)
%
% C.rpy(r,p,y) set camera rpy angles
% C.rpy(rpy)
%
% uv = C.project(P) return image coordinates for world points P
% uv = C.project(P, T) return image coordinates for world points P
% transformed by T prior to projection
%
% P is a list of 3D world points and the corresponding image plane points are
% returned in uv. Each point is represented by a column in P and uv.
%
% If P has 3 columns it is treated as a number of 3D points in world
% coordinates, one point per row.
%
% If POINTS has 6 columns, each row is treated as the start and end 3D
% coordinate for a line segment, in world coordinates.
%
% The optional arguments, T, represents a transformation that can be applied
% to the object data, P, prior to 'imaging'. The camera pose, C.Tcam, is also
% taken into account.
%
% uv = C.plot(P) display image coordinates for world points P
% uv = C.plot(P, T) isplay image coordinates for world points P transformed by T
%
% Points are displayed as a round marker. Lines are displayed as line segments.
% Optionally returns image plane coordinates uv.
%
% C.show
% C.show(name)
%
% Create a graphical camera with name, and pixel dimensions given by C.npix.
% Automatically called on first call to plot().
%
% SEE ALSO: Camera
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
% TODO:
% make a parent imaging class and subclass perspective, fisheye, panocam
% test for points in front of camera and set to NaN if not
% test for points off the image plane and set to NaN if not
% make clipping/test flags
classdef CatadioptricCamera < Camera
properties
k % radial distortion vector
model % projection model
maxangle % maximum elevation angle (above horizontal)
end
properties (SetAccess = private)
end
properties (GetAccess = private, SetAccess = private)
end
properties (GetAccess = private)
end
properties (Dependent = true, SetAccess = private)
end
methods
%
% Return a camera intrinsic parameter structure:
% focal length 8mm
% pixel size 10um square
% image size 1024 x 1024
% principal point (512, 512)
function c = CatadioptricCamera(varargin)
% invoke the superclass constructor
c = c@Camera(varargin{:});
c.type = 'FishEye';
if nargin == 0,
% default values
c.type = 'catadioptric';
c.k = 1;
c.model = 'equiangular';
c.name = 'catadioptric-default';
else
if isempty(c.k)
% compute k if not specified, so that hemisphere fits into
% image plane
r = min([(c.npix-c.pp).*c.rho, c.pp.*c.rho]);
switch c.model
case 'equiangular'
c.k = r / (pi/2 + c.maxangle);
case 'sine'
c.k = r;
case 'equisolid'
c.k = r / sin(pi/4);
case 'stereographic'
c.k = r / tan(pi/4);
r = c.k * tan(theta/2);
otherwise
error('unknown fisheye projection model');
end
end
end
end
function n = paramSet(c, args)
switch lower(args{1})
case 'maxangle'
c.maxangle = args{2}; n = 1;
case 'k'
c.k = args{2}; n = 1;
case 'projection'
c.model = args{2}; n = 1;
case 'default'
c.s = [10e-6, 10e-6]; % square pixels 10um side
c.npix = [1024, 1024]; % 1Mpix image plane
c.pp = [512, 512]; % principal point in the middle
c.limits = [0 1024 0 1024];
c.name = 'default';
r = min([(c.npix-c.pp).*c.s, c.pp.*c.s]);
c.k = 2*r/pi;
n = 0;
otherwise
error( sprintf('unknown option <%s>', args{count}));
end
end
function s = char(c)
s = sprintf('name: %s [%s]', c.name, c.type);
s = strvcat(s, sprintf( ' model: %s', c.model));
s = strvcat(s, sprintf( ' k: %-11.4g', c.k));
s = strvcat(s, char@Camera(c) );
end
% return field-of-view angle for x and y direction (rad)
function th = fov(c)
th = 2*atan(c.npix/2.*c.s / c.f);
end
% do the fisheye projection
function uv = project(c, P, varargin)
np = numcols(P);
opt.Tobj = [];
opt.Tcam = [];
[opt,arglist] = tb_optparse(opt, varargin);
if isempty(opt.Tcam)
opt.Tcam = c.T;
end
if isempty(opt.Tobj)
opt.Tobj = eye(4,4);
end
% transform all the points to camera frame
X = homtrans(inv(opt.Tcam) * opt.Tobj, P); % project them
R = colnorm(X(1:3,:));
phi = atan2( X(2,:), X(1,:) );
theta = acos( X(3,:) ./ R );
switch c.model
case 'equiangular'
r = c.k * theta;
case 'sine'
r = c.k * sin(theta);
case 'equisolid'
r = c.k * sin(theta/2);
case 'stereographic'
r = c.k * tan(theta/2);
otherwise
error('unknown fisheye projection model');
end
x = r .* cos(phi);
y = r .* sin(phi);
uv = [x/c.rho(1)+c.pp(1); y/c.rho(2)+c.pp(2)];
if c.noise
% add Gaussian noise with specified standard deviation
uv = uv + diag(c.noise) * randn(size(uv));
end
end
end % methods
end % class |
FishEyeCamera | Fish eye camera class
%
% C = camera default camera, 1024x1024, f=8mm, 10um pixels, camera at
% origin, optical axis is z-axis, u||x, v||y).
% C = camera(f, s, pp, npix, name)
% C = camera(0) f=sx=sy=1, u0=v0=0: normalized coordinates
%
% This camera model assumes central projection, that is, the focal point
% is at z=0 and the image plane is at z=f. The image is not inverted.
%
% The camera coordinate system is:
%
% 0------------> u X
% |
% |
% | + (principal point)
% |
% |
% v Y
% Z-axis is into the page.
%
% Object properties (read/write)
%
% C.f intrinsic: focal length
% C.s intrinsic: pixel size 2x1
% C.pp intrinsic: principal point 2x1
% C.np size of the virtual image plane (pixels) 2x1
%
% C.Tcam extrinsic: pose of the camera
% C.name name of the camera, used for graphical display
%
% Object properties (read only)
%
% C.su pixel width
% C.sv pixel height
% C.u0 principal point, u coordinate
% C.v0 principal point, v coordinate
%
% Object methods
%
% C.fov return camera half field-of-view angles (2x1 rads)
% C.K return camera intrinsic matrix (3x3)
% C.P return camera project matrix for camera pose (3x4)
% C.P(T) return camera intrinsic matrix for specified camera pose (3x4)
%
% C.rpy(r,p,y) set camera rpy angles
% C.rpy(rpy)
%
% uv = C.project(P) return image coordinates for world points P
% uv = C.project(P, T) return image coordinates for world points P
% transformed by T prior to projection
%
% P is a list of 3D world points and the corresponding image plane points are
% returned in uv. Each point is represented by a column in P and uv.
%
% If P has 3 columns it is treated as a number of 3D points in world
% coordinates, one point per row.
%
% If POINTS has 6 columns, each row is treated as the start and end 3D
% coordinate for a line segment, in world coordinates.
%
% The optional arguments, T, represents a transformation that can be applied
% to the object data, P, prior to 'imaging'. The camera pose, C.Tcam, is also
% taken into account.
%
% uv = C.plot(P) display image coordinates for world points P
% uv = C.plot(P, T) isplay image coordinates for world points P transformed by T
%
% Points are displayed as a round marker. Lines are displayed as line segments.
% Optionally returns image plane coordinates uv.
%
% C.show
% C.show(name)
%
% Create a graphical camera with name, and pixel dimensions given by C.npix.
% Automatically called on first call to plot().
%
% SEE ALSO: Camera
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
% TODO:
% make a parent imaging class and subclass perspective, fisheye, panocam
% test for points in front of camera and set to NaN if not
% test for points off the image plane and set to NaN if not
% make clipping/test flags
classdef FishEyeCamera < Camera
properties
k % radial distortion vector
model % projection model
end
properties (SetAccess = private)
end
properties (GetAccess = private, SetAccess = private)
end
properties (GetAccess = private)
end
properties (Dependent = true, SetAccess = private)
end
methods
%
% Return a camera intrinsic parameter structure:
% focal length 8mm
% pixel size 10um square
% image size 1024 x 1024
% principal point (512, 512)
function c = FishEyeCamera(varargin)
% invoke the superclass constructor
c = c@Camera(varargin{:});
c.type = 'FishEye';
if nargin == 0,
% default values
c.type = 'fisheye';
c.k = 1;
c.model = 'equiangular';
c.name = 'fisheye-default';
else
if isempty(c.k)
% compute k if not specified, so that hemisphere fits into
% image plane
r = min([(c.npix-c.pp).*c.rho, c.pp.*c.rho]);
switch c.model
case 'equiangular'
c.k = r / (pi/2);
case 'sine'
c.k = r;
case 'equisolid'
c.k = r / sin(pi/4);
case 'stereographic'
c.k = r / tan(pi/4);
r = c.k * tan(theta/2);
otherwise
error('unknown fisheye projection model');
end
end
end
end
function n = paramSet(c, args)
switch lower(args{1})
case 'k'
c.k = args{2}; n = 1;
case 'projection'
c.model = args{2}; n = 1;
case 'default'
c.rho = [10e-6, 10e-6]; % square pixels 10um side
c.npix = [1024, 1024]; % 1Mpix image plane
c.pp = [512, 512]; % principal point in the middle
c.limits = [0 1024 0 1024];
c.name = 'default';
r = min([(c.npix-c.pp).*c.rho, c.pp.*c.rho]);
c.k = 2*r/pi;
n = 0;
otherwise
error( sprintf('unknown option <%s>', args{count}));
end
end
function s = char(c)
s = sprintf('name: %s [%s]', c.name, c.type);
s = strvcat(s, sprintf( ' model: %s', c.model));
s = strvcat(s, sprintf( ' k: %-11.4g', c.k));
s = strvcat(s, char@Camera(c) );
end
% return field-of-view angle for x and y direction (rad)
function th = fov(c)
th = 2*atan(c.npix/2.*c.s / c.f);
end
% do the fisheye projection
function uv = project(c, P, varargin)
np = numcols(P);
opt.Tobj = [];
opt.Tcam = [];
[opt,arglist] = tb_optparse(opt, varargin);
if isempty(opt.Tcam)
opt.Tcam = c.T;
end
if isempty(opt.Tobj)
opt.Tobj = eye(4,4);
end
% transform all the points to camera frame
X = homtrans(inv(opt.Tcam) * opt.Tobj, P); % project them
R = colnorm(X(1:3,:));
phi = atan2( X(2,:), X(1,:) );
theta = acos( X(3,:) ./ R );
switch c.model
case 'equiangular'
r = c.k * theta;
case 'sine'
r = c.k * sin(theta);
case 'equisolid'
r = c.k * sin(theta/2);
case 'stereographic'
r = c.k * tan(theta/2);
otherwise
error('unknown fisheye projection model');
end
x = r .* cos(phi);
y = r .* sin(phi);
uv = [x/c.rho(1)+c.pp(1); y/c.rho(2)+c.pp(2)];
if c.noise
% add Gaussian noise with specified standard deviation
uv = uv + diag(c.noise) * randn(size(uv));
end
end
end % methods
end % class |
camcald | CAMCALD Compute camera calibration from data points
%
% C = CAMCALD(D) returns the camera calibration using a least squares
% technique. Input is a table of data points D with each row of the
% form [X Y Z U V] where (X,Y,Z) is the coordinate of a world point
% and [U,V] is the image plane coordinate of the corresponding point.
% C is a 3x4 camera calibration matrix.
%
% [C,E] = CAMCALD(D) as above but E is the maximum residual error after
% back substitution (in units of pixel).
%
% Note:
% - this method cannot handle lense distortion.
%
% See also CAMCALP, CAMERA, CAMCALT, INVCAMCAL.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function [C, resid] = camcald(XYZ, uv)
if numcols(XYZ) ~= numcols(uv)
error('must have same number of world and image-plane points');
end
n = numcols(XYZ);
if n < 6,
error('At least 6 points required for calibration');
end
%
% build the matrix as per Ballard and Brown p.482
%
% the row pair are one row at this point
%
A = [ XYZ' ones(n,1) zeros(n,4) -repmat(uv(1,:)', 1,3).*XYZ' ...
zeros(n,4) XYZ' ones(n,1) -repmat(uv(2,:)', 1,3).*XYZ' ];
%
% reshape the matrix, so that the rows interleave
%
A = reshape(A',11, n*2)';
if rank(A) < 11,
error('Rank deficient, perhaps points are coplanar or collinear');
end
B = reshape( uv, 1, n*2)';
C = A\B; % least squares solution
resid = max(max(abs(A * C - B)));
if resid > 1,
warning('Residual greater than 1 pixel');
end
fprintf('maxm residual %f pixels.\n', resid);
C = reshape([C;1]',4,3)'; |
FeatureMatch | % m display summary info about matches
% m.length number of matches
% m.inliers coordinates of inliers: u1 v1 u2 v2
% m.outliers coordinates of outliers: u1 v1 u2 v2
% m.plot display two source images side by side and overlay matches
% m.plot(''all'') " " show all matches (default)
% m.plot(''in'') " " only show inliers
% m.plot(''out'') " " only show outliers
% m.plot(''only'', n) " " only show n matches
% m.plot(''first'', n) " " only show first n matches
% m.plot(ls) " " with specified line style
% eg. m.plot(''out'', ''only'', 50, ''r'') show only 50 outliers as red lines
classdef FeatureMatch < handle
properties
xy_ % x1 y1 x2 y2 of corresponding points
distance_ % strength of match
inlier_ % NaN indeterminate
% true inlier
% false outlier
end
methods
function m = FeatureMatch(f1, f2, s)
if nargin == 0
return;
end
m.xy_ = [f1.u_ f1.v_ f2.u_ f2.v_]';
m.distance_ = s;
m.inlier_ = NaN;
end
function v = inlier(m)
v = m([m.inlier_] == true);
end
function v = outlier(m)
v = m([m.inlier_] == false);
end
function v = distance(m)
v = [m.distance_];
end
function display(m)
disp(' ');
disp([inputname(1), ' = '])
disp(' ');
if length(m) > 20
fprintf('%d corresponding points (listing suppressed)\n', length(m));
else
disp( char(m) );
end
end % display()
function s = char(matches)
s = '';
for m=matches
ss = sprintf('(%g, %g) <-> (%g, %g), dist=%f', ...
m.xy_, m.distance_);
switch m.inlier_
case true
ss = [ss ' +'];
case false
ss = [ss ' -'];
end
s = strvcat(s, ss);
end
end
function s = show(m)
s = sprintf('%d corresponding points\n', length(m));
in = [m.inlier_];
s = [s sprintf('%d inliers (%.1f%%)\n', ...
sum(in==true), sum(in==true)/length(m)*100)];
s = [s sprintf('%d outliers (%.1f%%)\n', ...
sum(in==false), sum(in==false)/length(m)*100) ];
end
function v = subset(m, n)
i = round(linspace(1, length(m), n));
v = m(i);
end
function s = p1(m, k)
xy = [m.xy_];
s = xy(1:2,:);
end
function s = p2(m, k)
xy = [m.xy_];
s = xy(3:4,:);
end
function s = p(m, k)
s = [m.xy_];
end
function plot(m, varargin)
try
ud = get(gca, 'UserData');
u0 = ud.u0;
catch
error('Current image is not a pair displayed by idisp');
end
w = u0(2);
xy = [m.xy_];
hold on
for k=1:numcols(xy),
plot([xy(1,k) xy(3,k)+w], xy([2 4],k), varargin{:});
end
hold off
figure(gcf);
end % plot
function [MM,rr] = ransac(m, func, varargin)
[M,in,resid] = ransac(func, [m.xy_], varargin{:});
% mark all as outliers
for i=1:length(m)
m(i).inlier_ = false;
end
for i=in
m(i).inlier_ = true;
end
if nargout >= 1
MM = M;
end
if nargout >= 2
rr = resid;
end
end
end
end |
ScalePointFeature | Point feature object
%
% A superclass for image corner features.
classdef ScalePointFeature < PointFeature
properties
scale_
end % properties
methods
function f = ScalePointFeature(varargin)
f = f@PointFeature(varargin{:}); % invoke the superclass constructor
end
function val = scale(features)
val = [features.scale_];
end
% accepts all the same options as imarker, first option must be the fill color
function plot_scale(features, varargin)
opt.display = {'circle', 'disk'};
[opt,arglist] = tb_optparse(opt, varargin);
holdon = ishold;
hold on
s = 1;
switch (opt.display)
case 'circle'
plot_circle([ [features.u_]; [features.v_] ], s*[features.scale_]', arglist{:});
case 'disk'
plot_circle([ [features.u_]; [features.v_] ], s*[features.scale_]', ...
'fillcolor', 'g', 'alpha', 0.2);
end
if ~holdon
hold off
end
end % plot
end % methods
end % classdef |
SurfPointFeature | SIFT Corner feature object
%
% A superclass for image corner features.
classdef SurfPointFeature < ScalePointFeature
properties
theta_
image_id_
end % properties
methods
function f = SurfPointFeature(varargin)
f = f@ScalePointFeature(varargin{:}); % invoke the superclass constructor
end
function val = theta_v(features)
val = [features.theta];
end
function val = theta(features)
val = [features.theta_];
end
function val = image_id(features)
val = [features.image_id_];
end
% accepts all the same options as imarker, first option must be the fill color
function plot_scale(features, varargin)
arglist = {};
argc = 1;
opt.display = 'circle';
while argc <= length(varargin)
switch lower(varargin{argc})
case 'circle'
opt.display = varargin{argc};
case 'clock'
opt.display = varargin{argc};
case 'disk'
opt.display = varargin{argc};
case 'arrow'
opt.display = varargin{argc};
otherwise
arglist = [arglist varargin(argc)];
end
argc = argc + 1;
end
holdon = ishold;
hold on
s = 20/sqrt(pi); % circle of same area as 20s x 20s square support region
switch (opt.display)
case 'circle'
plot_circle([ [features.u_]; [features.v_] ], s*[features.scale_]', arglist{:});
case 'clock'
plot_circle([ [features.u_]; [features.v_] ], s*[features.scale_]', arglist{:});
for f=features
plot([f.u_, f.u_+s*f.scale_*cos(f.theta_)], ...
[f.v_, f.v_+s*f.scale_*sin(f.theta_)], ...
arglist{:});
end
case 'disk'
plot_circle([ [features.u_]; [features.v_] ], s*[features.scale_]', ...
'fillcolor', 'g', 'alpha', 0.2);
case 'arrow'
for f=features
quiver(f.u_, f.v_, s*f.scale_.*cos(f.theta_), ...
s*f.scale_.*sin(f.theta_), arglist{:});
end
end
if ~holdon
hold off
end
end % plot
function [m,corresp] = match(f1, f2, varargin)
opt.thresh = 0.05;
opt = tb_optparse(opt, varargin);
% Put the landmark descriptors in a matrix
%D1 = reshape([Ipts1.descriptor],64,[]);
%D2 = reshape([Ipts2.descriptor],64,[]);
D1 = f1.descriptor;
D2 = f2.descriptor;
% Find the best matches
err=zeros(1,length(f1));
cor1=1:length(f1);
cor2=zeros(1,length(f1));
for i=1:length(f1),
distance=sum((D2-repmat(D1(:,i),[1 length(f2)])).^2,1);
[err(i),cor2(i)]=min(distance);
end
% Sort matches on vector distance
[err, ind]=sort(err);
cor1=cor1(ind);
cor2=cor2(ind);
m = [];
cor = [];
for i=1:length(f1)
k1 = cor1(i);
k2 = cor2(i);
mm = FeatureMatch(f1(k1), f2(k2), err(i));
m = [m mm];
cor(:,i) = [k1 k2]';
end
if isnumeric(opt.thresh)
thresh = opt.thresh;
elseif strcmp('median', opt.thresh)
thresh = median(err);
end
if ~isempty(opt.thresh)
k = err > thresh;
cor(:,k) = [];
m(k) = [];
end
if nargout > 1
corresp = cor;
end
%
% [matches,dist,dist2] = closest([f1.descriptor], [f2.descriptor]);
% matches = [1:length(f1); matches];
%
% % delete matches where distance of closest match is greater than
% % 0.7 of second closest match
% k = dist > 0.7 * dist2;
%
% matches(:,k) = [];
% dist(k) = [];
%
% % dist is a 1xM matrix of distance between the matched features, low is good.
%
% % matches is a 2xM matrix, one column per match, each column
% % is the index of the matching feature in images 1 and 2
%
% % sort into increasing distance
% [z,k] = sort(dist, 'ascend');
% matches = matches(:,k);
% dist = dist(:,k);
%
% m = [];
% cor = [];
%
% for i=1:numcols(matches),
% k1 = matches(1,i);
% k2 = matches(2,i);
% mm = FeatureMatch(f1(k1), f2(k2), dist(i));
% m = [m mm];
% cor(:,i) = [k1 k2]';
% end
%
% if nargout > 1
% corresp = cor;
% end
end
end % methods
methods(Static)
% the MEX functions live in a private subdirectory, so these static methods
% provide convenient access to them
function Ipts = surf(varargin)
Ipts = OpenSurf(varargin{:});
end
end
end % classdef |
PointFeature | PointCorner feature object
%
% A superclass for image corner features.
classdef PointFeature < handle
properties %(GetAccess=protected)%, Hidden=true)
u_ % feature x-coordinates
v_ % feature y-coordinates
strength_
descriptor_
end % properties
methods
function f = PointFeature(u, v, strength)
if nargin == 0
return;
end
if nargin >= 2
f.u_ = u;
f.v_ = v;
end
if nargin == 3
f.strength_ = strength;
end
end
function val = u(f)
val = [f.u_];
end
function val = v(f)
val = [f.v_];
end
function val = p(f)
val = [[f.u_]; [f.v_]];
end
function val = strength(f)
val = [f.strength_];
end
function val = descriptor(f)
val = [f.descriptor_];
end
function display(f)
disp(' ');
disp([inputname(1), ' = '])
disp(' ');
if length(f) > 20
fprintf('%d features (listing suppressed)\n Properties:', length(f));
for property=fieldnames(f)'
fprintf(' %s', property{1}(1:end-1));
end
fprintf('\n');
else
disp( char(f) );
end
end % display()
function ss = char(features)
ss = [];
for i=1:length(features)
f = features(i);
% display the coordinate
s = sprintf(' (%g,%g)', f.u_, f.v_);
% display the other properties
for property=fieldnames(f)'
prop = property{1}; % convert from cell array
switch prop
case {'u_', 'v_', 'descriptor_'}
continue;
otherwise
val = getfield(f, prop);
if ~isempty(val)
s = strcat(s, [sprintf(', %s=', prop(1:end-1)), num2str(val, ' %g')]);
end
end
end
% do the descriptor last
val = getfield(f, 'descriptor_');
if ~isempty(val)
if length(val) == 1
% only list scalars or shortish vectors
s = strcat(s, [', descrip=', num2str(val, ' %g')]);
elseif length(val) < 4
% only list scalars or shortish vectors
s = strcat(s, [', descrip=(', num2str(val', ' %g'), ')']);
else
s = strcat(s, ', descrip= ..');
end
end
ss = strvcat(ss, s);
end
end
function val = uv(features)
val = [[features.u]; [features.v]];
end
% f.plot()
% f.plot(linespec)
function plot(features, varargin)
holdon = ishold;
hold on
if nargin == 1
varargin = {'ws'};
end
for i=1:length(features)
plot(features(i).u_, features(i).v_, varargin{:});
end
if ~holdon
hold off
end
end % plot
function s = distance(f1, f2)
for i=1:length(f2)
s(i) = norm(f1.descriptor-f2(i).descriptor);
end
end
function s = ncc(f1, f2)
for i=1:length(f2)
s(i) = dot(f1.descriptor,f2(i).descriptor);
end
end
function [m,corresp] = match(f1, f2)
% TODO: extra args for distance measure, could be ncc, pass through to closest
% allow threshold, percentage of max
[corresp, dist] = closest([f1.xy_], [f2.xy_]);
% sort into increasing distance
[z,k] = sort(dist, 'ascend');
corresp = corresp(:,k);
dist = dist(:,k);
m = [];
cor = [];
for i=1:numcols(corresp),
k1 = i;
k2 = corresp(i);
mm = FeatureMatch(f1(k1), f2(k2), dist(i));
m = [m mm];
cor(:,i) = [k1 k2]';
end
if nargout > 1
corresp = cor;
end
end
end % methods
end % classdef |
iscalemax | Find maxima in scale space image sequence
%
% F = ISCALEMAX(L, S) returns a vector of ScalePointFeature objects from the
% scale space image sequence L which is NxMxD. S is a vector of scale values
% corresponding to each plane of L
%
% Notes::
% - feautures are sorted into descending feature strength
%
% See also ISCALESPACE, ScalePointFeature.
function features = iscalemax(L, sscale)
allmax = zeros(size(L));
nbmax = zeros(size(L));
se_all = ones(3,3);
se_nb = se_all; se_nb(2,2) = 0;
% get maxima at each level
for k=1:size(L,3)
allmax(:,:,k) = imorph(abs(L(:,:,k)), se_all, 'max', 'replicate');
nbmax(:,:,k) = imorph(abs(L(:,:,k)), se_nb, 'max', 'replicate');
end
z = zeros(size(L,1), size(L,2));
fc = 1;
strength = zeros(size(L,1), size(L,2));
for k=2:size(L,3)-1 % maxima cant be at either end of the scale range
s = abs(L(:,:,k));
corners = find( s > nbmax(:,:,k) & s > allmax(:,:,k-1) & s > allmax(:,:,k+1) );
for corner=corners'
[y,x] = ind2sub(size(s), corner);
features(fc) = ScalePointFeature(x, y, s(corner));
features(fc).scale_ = sscale(k)*sqrt(2);
fc = fc+1;
end
end
% sort into descending order of strength
[z,k] = sort(-features.strength);
features = features(k); |
iscalespace | Return scale-space image sequence
%
% [G,L,S] = ISCALESPACE(IM, D, SIGMA) returns a scale space image sequence of
% length D derived from IM (HxW). The standard deviation of the smoothing
% Gaussian is SIGMA. At each scale step the variance of the Gaussian increases
% by SIGMA^2. The first step in the sequence is the original image.
%
% L (HxWxN) is the absolute value of the Laplacian of the scale sequence,
% G (HxWxN) is the scale sequence, and S (Nx1) is the vector of scales
% corresponding to each step of the sequence.
%
% [G,L,S] = ISCALESPACE(IM, N) as above but SIGMA=1.
%
% Notes::
% - The Laplacian is computed by the difference of adjacent Gaussians.
%
% See also ISCALEMAX, ISMOOTH, ILAPLACE.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function [G, L, S] = scalespace(im, n, sigma)
if nargin < 3
sigma = 1;
end
g(:,:,1) = im;
scale = 0.5;
scales = scale;
for i=1:n-1,
im = ismooth(im, sigma);
scale = sqrt(scale^2 + sigma^2);
scales = [scales; scale];
g(:,:,i+1) = im;
lap(:,:,i) = scale^2 * ( g(:,:,i+1) - g(:,:,i) );
end
% return results as requested
if nargout > 2,
S = scales;
end
if nargout > 0,
G = g;
end
if nargout > 1,
L = lap;
end |
isurf | SURF feature extractor
%
% SF = ISURF(IM, OPTIONS) returns a vector of Surf objects
% representing scale and rotationally invariant interest points in the
% image IM.
%
% The SurfPointFeature object has many properties including:
% u horizontal coordinate
% v vertical coordinate
% strength feature strength
% descriptor feature descriptor (64x1 or 128x1)
% sigma feature scale
% theta feature orientation (rad)
%
% Options::
% 'nfeat',N set the number of features to return (default Inf)
% 'extended' return 128-element descriptor (default 64)
% 'upright' dont compute rotation invariance
% 'suppress',R set the suppression radius (default 0)
%
% Notes::
% - features are returned in descending strength order
% - wraps a function by ??
% - if IM is NxMxP it is taken as an image sequence and F is a cell array whose
% elements are feature vectors for the corresponding image in the sequence.
%
% Reference:
% "SURF: Speeded Up Robust Features", Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool,
% Computer Vision and Image Understanding (CVIU), Vol. 110, No. 3, pp. 346--359, 2008
%
% See also SurfPointFeature, ISIFT, ICORNER.
function features = isurf(im, varargin)
opt.suppress = 0;
opt.nfeat = Inf;
opt.extended = true;
opt.upright = false;
[opt,arglist] = tb_optparse(opt, varargin);
if iscell(im)
% images provided as a cell array, return a cell array
% of SIFT object vectors
fprintf('extracting SIFT features for %d greyscale images\n', length(im));
features = {};
for i=1:length(im)
sf = isurf(im{i}, 'setopt', opt);
for j=1:length(sf)
sf(j).image_id = i;
end
features{i} = sf;
fprintf('.');
end
fprintf('\n');
return
end
% convert color image to greyscale
if ndims(im) ==3 && size(im, 3) == 3
im = imono(im);
end
if ndims(im) > 2
% TODO sequence of color images..
% images provided as an array, return a cell array
% of SURF object vectors
if opt.verbose
fprintf('extracting SURF features for %d greyscale images\n', size(im,3));
end
features = {};
for i=1:size(im,3)
sf = isurf(im(:,:,i), 'setopt', opt);
for j=1:length(sf)
sf(j).image_id_ = i;
end
features{i} = sf;
fprintf('.');
end
if opt.verbose
fprintf('\n');
end
return
end
% do SURF using a static method that wraps the implementation from:
% OpenSURF for Matlab
% written by D.Kroon University of Twente (July 2010)
% based on the C++ implementation by Chris Evans
%options.extended = true;
options.extended = false;
options.tresh = 0.0002;
Ipts = SurfPointFeature.surf(im, options, arglist{:});
% Ipts is a structure array with elements x, y, scale, orientation, descriptor
fprintf('%d corners found (%.1f%%), ', length(Ipts), ...
length(Ipts)/prod(size(im))*100);
% sort into descending order of corner strength
[z,k] = sort([Ipts.strength], 'descend');
Ipts = Ipts(:,k);
% allocate storage for the objects
n = min(opt.nfeat, length(Ipts));
features = [];
i = 1;
while i<=n
if i > length(Ipts)
break;
end
% enforce separation between corners
% TODO: strategy of Brown etal. only keep if 10% greater than all within radius
if (opt.suppress > 0) && (i>1)
d = sqrt( ([features.v]'-Ipts(i).y).^2 + ([features.u]'-Ipts(i).x).^2 );
if min(d) < opt.suppress
continue;
end
end
f = SurfPointFeature(Ipts(i).x, Ipts(i).y, Ipts(i).strength);
f.scale_ = Ipts(i).scale;
f.theta_ = Ipts(i).orientation;
f.descriptor_ = cast(Ipts(i).descriptor, 'single');
features = [features f];
i = i+1;
end
fprintf(' %d corner features saved\n', i-1); |
icorner | Classical corner detector
%
% F = ICORNER(IM, OPTIONS) returns a vector of PointFeature objects describing
% the detected corner features. This is a non-scale space detector and by
% default the Harris method is used. If IM is an image sequence a cell array
% of feature vectors is returned.
%
% The PointFeature object has many properties including:
% u horizontal coordinate
% v vertical coordinate
% strength corner strength
% descriptor corner descriptor (vector)
%
% Options::
% 'cmin',CM minimum corner strength
% 'cminthresh',CT minimum corner strenght as a fraction of maximum corner strength
% 'edgegap',EG dont return features closer than EG to the edge of image
% 'suppress',R dont return a feature closer than R pixels to an earlier feature
% 'nfeat',N return the N strongest corners (default Inf)
% 'detector',D choose the detector where D is 'harris' (default), 'noble' or 'klt'
% 'sigma',S kernel width for smoothing
% 'deviv',D kernel for gradient (default kdgauss(2))
% 'k',K set the value of k for Harris detector
% 'patch',P use a PxP patch as the feature vector
% 'color' specify that IM is a color image not a sequence
%
% Notes::
% - corners are processed in order from strongest to weakest.
% - the function stops when:
% - the corner strength drops below .cmin
% - the corner strenght drops below P.cMinThresh * strongest corner
% - the list of corners is exhausted
% - features are returned in descending strength order
% - if IM has more than 2 dimensions it is either a color image or a sequence
% - if IM is NxMxP it is taken as an image sequence and F is a cell array whose
% elements are feature vectors for the corresponding image in the sequence.
% - if IM is NxMx3 it is taken as a sequence unless the option 'color' is used
% - if IM is NxMx3xP it is taken as a sequence of color images and F is a cell
% array whose elements are feature vectors for the corresponding color image
% in the sequence.
% - the default descriptor is a vector [Ix Iy Ixy] which are the unique
% elements of the structure tensor.
%
% Reference:
% "A combined corner and edge detector", C.G. Harris and M.J. Stephens,
% Proc. Fourth Alvey Vision Conf., Manchester, pp 147-151, 1988.
%
% See also PointFeature, ISURF.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function [features, corner_strength] = icorner(im, varargin)
% TODO, can handle image sequence, return 3D array of corner_strength if requested
% and cell array of corner vectors
% parse options into parameter struct
opt.k = 0.04;
opt.deriv = kdgauss(2);
opt.cmin = 0;
opt.cminthresh = 0.0;
opt.edgegap = 2;
opt.nfeat = Inf;
opt.sigma_i = 2;
opt.patch = 0;
opt.detector = {'harris', 'noble', 'st'};
opt.color = false;
opt.suppress = 0;
opt.nfeat = 100;
[opt,arglist] = tb_optparse(opt, varargin);
if opt.patch > 0
opt.edgegap = opt.patch;
end
if iscell(im)
% images provided as a cell array, return a cell array
% of corner object vectors
if opt.verbose
fprintf('extracting corner features for %d images\n', length(im));
end
features = {};
for i=1:length(im)
f = icorner(im{i}, 'setopt', opt);
for j=1:length(f)
f(j).image_id = i;
end
features{i} = f;
fprintf('.');
end
if opt.verbose
fprintf('\n');
end
return
end
if ndims(im) > 2
% images provided as an array, return a cell array
% of corner object vectors
if ndims(im) == 4 && size(im,3) == 3 && opt.color
fprintf('extracting corner features for color %d images\n', size(im,4));
features = {};
% sequence of color images
for i=1:size(im,k)
f = icorner(im(:,:,:,i), 'setopt', opt);
for j=1:length(f)
%f(j).image_id_ = i;
end
features{i} = f;
fprintf('.');
end
fprintf('\n');
return
elseif ndims(im) == 3 && ~opt.color
fprintf('extracting corner features for grey %d images\n', size(im,3));
features = {};
% sequence of grey images
for i=1:size(im,3)
f = icorner(im(:,:,i), 'setopt', opt);
for j=1:length(f)
%f(j).image_id_ = i;
end
features{i} = f;
fprintf('.');
end
fprintf('\n');
return
end
end
if ndims(im) == 3 & opt.color
R = double(im(:,:,1));
G = double(im(:,:,2));
B = double(im(:,:,3));
Rx = conv2(R, opt.deriv, 'same');
Ry = conv2(R, opt.deriv', 'same');
Gx = conv2(G, opt.deriv, 'same');
Gy = conv2(G, opt.deriv', 'same');
Bx = conv2(B, opt.deriv, 'same');
By = conv2(B, opt.deriv', 'same');
Ix = Rx.^2+Gx.^2+Bx.^2;
Iy = Ry.^2+Gy.^2+By.^2;
Ixy = Rx.*Ry+Gx.*Gy+Bx.*By;
else
% compute horizontal and vertical gradients
im = double(im);
ix = conv2(im, opt.deriv, 'same');
iy = conv2(im, opt.deriv', 'same');
Ix = ix.*ix;
Iy = iy.*iy;
Ixy = ix.*iy;
end
% smooth them
if opt.sigma_i > 0
Ix = ismooth(Ix, opt.sigma_i);
Iy = ismooth(Iy, opt.sigma_i);
Ixy = ismooth(Ixy, opt.sigma_i);
end
[nr,nc] = size(Ix);
npix = nr*nc;
% computer cornerness
switch opt.detector
case 'harris'
cornerness = (Ix .* Iy - Ixy.^2) - opt.k * (Ix + Iy).^2;
case 'noble'
cornerness = (Ix .* Iy - Ixy.^2) ./ (Ix + Iy);
case 'st'
cornerness = zeros(size(Ix));
for i=1:npix
lambda = eig([Ix(i) Ixy(i); Ixy(i) Iy(i)]);
cornerness(i) = min(lambda);
end
end
% compute maximum value around each pixel
cmax = imorph(cornerness, [1 1 1;1 0 1;1 1 1], 'max');
% if pixel exceeds this, its a local maxima, find index
cindex = find(cornerness > cmax);
fprintf('%d corners found (%.1f%%), ', length(cindex), ...
length(cindex)/npix*100);
% remove those near edges
[y, x] = ind2sub(size(cornerness), cindex);
e = opt.edgegap;
k = (x>e) & (y>e) & (x < (nc-e)) & (y < (nr-e));
cindex = cindex(k);
% corner strength must exceed an absolute minimum
k = cornerness(cindex) < opt.cmin;
cindex(k) = [];
% sort into descending order
cval = cornerness(cindex); % extract corner values
[z,k] = sort(cval, 'descend'); % sort into descending order
cindex = cindex(k)';
cmax = cornerness( cindex(1) ); % take the strongest feature value
% corner strength must exceed a fraction of the maximum value
k = cornerness(cindex)/cmax < opt.cminthresh;
cindex(k) = [];
% allocate storage for the objects
n = min(opt.nfeat, numcols(cindex));
features = [];
i = 1;
while i <= n
if i > length(cindex)
break;
end
K = cindex(i);
c = cornerness(K);
% get the coordinate
[y, x] = ind2sub(size(cornerness), K);
% enforce separation between corners
% TODO: strategy of Brown etal. only keep if 10% greater than all within radius
if (opt.suppress > 0) && (i>1)
d = sqrt( sum((features.v'-y).^2 + (features.u'-x).^2) );
if min(d) < opt.suppress
continue;
end
end
% ok, this one is for keeping
f = PointFeature(x, y, c);
if opt.patch == 0
f.descriptor_ = [Ix(K) Iy(K) Ixy(K)]';
else
% if opt.patch is finite, then return a vector which is the local image
% region as a vector, zero mean, and normalized by the norm.
% the dot product of this with another descriptor is the ZNCC similarity measure
w2 = opt.patch;
d = im(y-w2:y+w2,x-w2:x+w2);
d = d(:);
d = d - mean(d);
f.descriptor_ = cast( d / norm(d), 'single');
end
features = [features f];
i = i+1;
end
fprintf(' %d corner features saved\n', i-1);
% sort into descending order of strength
[z,k] = sort(-features.strength);
features = features(k);
if nargout > 1
corner_strength = cornerness;
end |
istereo | Stereo matching
%
% D = ISTEREO(IML, IMR, W, RANGE, OPTIONS) returns a disparity image the same
% size as the input images IML and IMR and the value at each pixel is the
% horizontal shift of the corresponding pixel in IML as observed in IMR.
%
% IML and IMR are the left- and right-images of a stereo pair, of either
% double or uint8 class.
%
% W is the size of the matching window, which can be a scalar for WxW or a
% 2-vector [WX WY] for a WXxWY window.
%
% RANGE is the disparity search range, which can be a scalar for disparities in
% the range 0 to RANGE, or a 2-vector [DMIN DMAX] for searches in the range
% DMIN to DMAX.
%
% [D,SIM] = ISTEREO(IML, IMR, W, RANGE, OPTIONS) as above but returns SIM
% which is the same size as D and the elements are the peak matching score
% for the corresponding elements of D. For the default matching metric ZNCC
% this varies between -1 (very bad) to +1 (perfect).
%
% [D,SIM,DSI] = ISTEREO(IML, IMR, W, RANGE, OPTIONS) as above but returns DSI
% which is the disparity space image. If IML and IMR are NxM and range
% specifies D % disparity values then DSI is NxMxD where the I'th plane is the
% similarity of IML to IMR shifted by ????
%
% [D,SIM,P] = ISTEREO(IML, IMR, W, RANGE, OPTIONS) if the 'interp' option is
% given then disparity is estimated to sub-pixel precision using quadratic
% interpolation. In this case D is the interpolated disparity and P is
% a structure with elements A, B, dx. The interpolation polynomial is
% s = Ad^2 + Bd + C where s is the similarity score and d is disparity relative
% to the integer disparity at which s is maximum. P.A and P.B are matrices the
% same size as D whose elements are the per pixel values of the interpolation
% polynomial coefficients. P.dx is the peak of the polynomial with respect
% to the integer disparity at which s is maximum (in the range -0.5 to +0.5).
%
% Options::
% 'metric',M string that specifies the similarity metric to use which is
% one of 'zncc' (default), 'ncc', 'ssd' or 'sad'.
% 'interp' enable subpixel interpolation and D contains non-integer
% values (default false)
%
% Notes::
% - For both output images the pixels within a half-window dimension of the edges
% will not be valid and are set to NaN.
% - SIM = max(DSI, 3)
%
% See also STVIEW.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function [disp,sim, o3] = istereo(L, R, drange, h, varargin)
% TODO: score cube is float, can return it
opt.metric = 'zncc';
opt.interp = false;
opt = tb_optparse(opt, varargin);
% ensure images are greyscale
L = imono(L);
R = imono(R);
if length(drange) > 2
vshift = drange(3);
if vshift > 0
L = L(vshift:end,:);
R = R(1:end-vshift,:);
else
vshift = -vshift;
L = L(1:end-vshift,:);
R = R(vshift:end,:);
end
end
% compute the score cube, 3rd dimension is disparity
DSI = stereo_match(L, R, 2*h+1, drange(1:2), opt.metric);
% best value along disparity dimension is the peak
% s best score
% d disparity at which it occurs
%
% both s and d are matrices same size as L and R.
if strcmp(opt.metric, 'sad') | strcmp(opt.metric, 'ssd')
[s,d] = min(DSI, [], 3);
else
[s,d] = max(DSI, [], 3);
end
d(isnan(s)) = NaN;
if opt.interp
% interpolated result required
% get number of pixels and disparity range
npix = prod(size(L));
if length(drange) == 1,
ndisp = drange + 1;
else
dmin = min(drange);
dmax = max(drange);
ndisp = dmax - dmin + 1;
end
% find all disparities that are not at either end of the range, we need
% a point on either side to interpolate them
valid = (d>1) & (d 0,
disp = d;
end
if nargout > 1,
sim = s;
end
if nargout > 2
if opt.interp
o3.A = A;
o3.B = B;
o3.dx = dx;
else
o3 = DSI;
end
end |
anaglyph | Convert stereo images to an anaglyph image
%
% A = ANAGLYPH(LEFT, RIGHT) is anaglyph image where the two images of
% stereo pair are coded are combined in a single image using two different
% colors. By default the left image is red, and the right image is cyan.
%
% A = ANAGLYPH(LEFT, RIGHT, COLOR) as above but the string COLOR describes
% the color coding as a string with 2 letters, the first for left, the second
% for right, and each is one of:
% r red
% g green
% b green
% c cyan
% m magenta
% ag = anaglyph(left, right, colors, disp)
% ag = anaglyph(stereopair, colors, disp)
%
% A = ANAGLYPH(LEFT, RIGHT, COLOR, DISP) as above but allows for disparity
% correction. If DISP is positive the disparity is increased, if negative it
% is reduced. This is achieved by trimming the images. Use this option to
% make the images more natural/comfortable to view, useful if the images were
% achieved with a non-human stereo baseline or field of view.
%
% See also STVIEW.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function anaglyph = anaglyph(left, right, colors, disp)
if nargin < 3,
colors = 'rc';
end
if nargin < 4,
disp = 0;
end
% ensure the images are greyscale
left = imono(left);
right = imono(right);
[height,width] = size(left);
if disp > 0,
left = left(:,1:width-disp);
right = right(:,disp+1:end);
end
if disp < 0,
disp = -disp;
left = left(:,disp+1:end);
right = right(:,1:width-disp);
end
ag = zeros([size(left) 3]);
ag = ag_insert(ag, left, colors(1));
ag = ag_insert(ag, right, colors(2));
if nargout > 0,
aglyph = ag;
else
if isa(left, 'uint8'),
ag = ag / 255;
end
image(ag);
end
function out = ag_insert(in, im, c)
out = in;
% map single letter color codes to image planes
switch c,
case 'r'
out(:,:,1) = im; % red
case 'g'
out(:,:,2) = im; % green
case 'b'
% blue
out(:,:,3) = im;
case 'c'
out(:,:,2) = im; % cyan
out(:,:,3) = im;
case 'm'
out(:,:,1) = im; % magenta
out(:,:,3) = im;
case 'o'
out(:,:,1) = im; % orange
out(:,:,2) = im;
end |
stdisp | Display stereo pair
%
% STDISP(L, R) displays the stereo image pair L and R in adjacent windows.
%
% Two cross-hairs are created. Clicking a point in the left image positions
% black cross hair at the same pixel coordinate in the right image. Clicking
% the corresponding world point in the right image sets the green crosshair
% and displays the disparity.
%
% See also ISTEREO.
function stdisp(L, R)
% display the images side by side
idisp([L R], 'nogui');
% initial cross hair location
Y = 100;
X = 100;
% create the 3 lines segments and stash in user data on the axis
ud.w = size(L,2); % width of left image
ud.hline = line('XData',get(gca,'XLim'),'YData',[Y Y], ...
'Tag','Horizontal Cursor');
ud.vline_l = line('XData',[X X],'YData',get(gca,'YLim'), ...
'Tag','Vertical Cursor');
ud.vline_r = line('XData',[X+ud.w X+ud.w],'YData',get(gca,'YLim'), ...
'Tag','Vertical Cursor');
ud.vline_r2 = line('XData',[X+ud.w X+ud.w],'YData',get(gca,'YLim'), ...
'Tag','Vertical Cursor', 'color', 'g');
ud.panel = uicontrol(gcf, ...
'style', 'text', ...
'units', 'norm', ...
'pos', [.5 .935 .48 .05], ...
'background', [1 1 1], ...
'HorizontalAlignment', 'left', ...
'string', ' Machine Vision Toolbox for Matlab ' ...
);
set(gca, 'UserData', ud);
% Set the WindowButtonFcn of the figure
set(gcf,'WindowButtonDownFcn', @buttonDown,...
'WindowButtonUpFcn',@buttonUp);
end
function moveCursor(src, event)
ud = get(gca, 'UserData');
cp = get(gca,'CurrentPoint');
% cp = [xfront yfront xfront; xback yback zback]
if cp(1,1) < ud.w
set(ud.hline, 'YData', [cp(1,2) cp(1,2)]);
set(ud.vline_l, 'XData', [cp(1,1) cp(1,1)]);
set(ud.vline_r, 'XData', ud.w+[cp(1,1) cp(1,1)]);
else
set(ud.vline_r2, 'XData', [cp(1,1) cp(1,1)]);
xl = get(ud.vline_l, 'XData');
%fprintf('d = %f\n', cp(1,1) - xl(1) - ud.w);
set(ud.panel, 'string', sprintf('d = %f\n', xl(1) + ud.w - cp(1,1)));
end
end
function buttonDown(src, event)
set(gcf, 'WindowButtonMotionFcn',@moveCursor);
moveCursor(src, event);
%disp('down');
end
function buttonUp(src, event)
%disp('up');
set(gcf, 'WindowButtonMotionFcn','');
end |
irectify | Rectify stereo image pair
%
% [OUT1,OUT2] = IRECTIFY(F, M, IM1, IM2) returns a rectified pair of images
% corresponding to IM1 and IM2. F is the fundamental matrix relating the two
% views and M is a Match object containing point correspondences between the
% images.
%
% [OUT1,OUT2,H1,H2] = IRECTIFY(F, M, IM1, IM2) as above but also returns
% the homographies H1 and H2 that map IM1 to OUT1 and IM2 to OUT2 respectively.
%
% See also ISTEREO, CentralCamera.
function [Img1_new, Img2_new, H12,H21] = irectify(F, m, Img1, Img2)
% http://se.cs.ait.ac.th/cvwiki/matlab:tutorial:rectification
F12 = F';
[rows,cols,depth] = size(Img1);
% Get homographies.
x1 = e2h( m.p1 );
x2 = e2h( m.p2 );
[H12,H21,bSwap] = rectify_homographies( F12, x1, x2, rows, cols );
[w1,off1] = homwarp(H12, Img1, 'full');
[w2,off2] = homwarp(H21, Img2, 'full');
% fix the vertical alignment of the images by padding
dy = off1(2) - off2(2);
if dy < 0
w1 = ipad(w1, 'b', -dy);
w2 = ipad(w2, 't', -dy);
else
w1 = ipad(w1, 't', dy);
w2 = ipad(w2, 'b', dy);
end
[w1,w2] = itrim(w1, w2);
if nargout == 0
stview(w1, w2)
else
Img1_new = w1;
Img2_new = w2;
end
%-----------------------------------------------------------------------------
function [H1,H2,bSwap] = rectify_homographies( F, x1, x2, rows, cols )
% F: a fundamental matrix
% x1 and x2: corresponding points such that x1_i' * F * x2_i = 0
% Initialize
H1 = [];
H2 = [];
bSwap = 0;
% Center of image
cy = round( rows/2 );
cx = round( cols/2 );
% Fix F to be rank 2 to numerical accuracy
[U,D,V] = svd( F );
D(3,3) = 0;
F = U*D*V';
% Get epipole. e12 is the epipole in image 1 for camera 2.
e12 = null( F' ); % Epipole in image 1 for camera 2
e21 = null( F ); % Epipole in image 2 for camera 1
% Put epipoles in front of camera
if e12 < 0, e12 = -e12; end;
if e21 < 0, e21 = -e21; end;
% Make sure the epipoles are inside the images
check_epipoles_in_image( e12, e21, rows, cols );
% Check that image 1 is to the left of image 2
% if e12(1)/e12(3) < cx
% fprintf( 1, 'Swapping left and right images...\n' );
% tmp = e12;
% e12 = e21;
% e21 = tmp;
% F = F';
% bSwap = 1;
% end;
% Now we have
% F' * e12 = 0,
% F * e21 = 0,
% Let's get the rectifying homography Hprime for image 1 first
Hprime = map_to_infinity( e12, cx, cy );
e12_new = Hprime * e12;
% Normalize Hprime so that Hprime*eprime = (1,0,0)'
Hprime = Hprime / e12_new(1);
e12_new = Hprime * e12;
fprintf( 1, 'Epipole 1/2 mapped to infinity: (%g, %g, %g)\n', e12_new );
% Get canonical camera matrices for F12 and compute H0, one possible
% rectification homography for image 2
[P,Pprime] = get_canonical_cameras( F );
M = Pprime(:,1:3);
H0 = Hprime * M;
% Test that F12 is a valid F for P,Pprime
test_p_f( P, Pprime, F );
% Now we need to find H so that the epipolar lines match
% each other, i.e., inv(H)' * l = inv(Hprime)' * lprime
% and the disparity is minimized, i.e.,
% min \sum_i d(H x_i, Hprime xprime_i)^2
% Transform data initially according to Hprime (img 1) and H0 (img 2)
x1hat = Hprime * x1;
x1hat = x1hat ./ repmat( x1hat(3,:), 3, 1 );
x2hat = H0 * x2;
x2hat = x2hat ./ repmat( x2hat(3,:), 3, 1 );
rmse_x = sqrt( mean( (x1hat(1,:) - x2hat(1,:) ).^2 ));
rmse_y = sqrt( mean( (x1hat(2,:) - x2hat(2,:) ).^2 ));
fprintf( 1, 'Before Ha, RMSE for corresponding points in Y: %g X: %g\n', ...
rmse_y, rmse_x );
% Estimate [ a b c ; 0 1 0 ; 0 0 1 ] aligning H, Hprime
n = size(x1,2);
A = [ x2hat(1,:)', x2hat(2,:)', ones(n,1) ];
b = x1hat(1,:)';
abc = A\b;
HA = [ abc' ; 0 1 0 ; 0 0 1 ];
H = HA*H0;
x2hat = H * x2;
x2hat = x2hat ./ repmat( x2hat(3,:), 3, 1 );
rmse_x = sqrt( mean(( x1hat(1,:) - x2hat(1,:) ).^2 ));
rmse_y = sqrt( mean(( x1hat(2,:) - x2hat(2,:) ).^2 ));
fprintf( 1, 'After Ha, RMSE for corresponding points in Y: %g X: %g\n', ...
rmse_y, rmse_x );
% Return the homographies as appropriate
if bSwap
H1 = H;
H2 = Hprime;
else
H1 = Hprime;
H2 = H;
end;
%-----------------------------------------------------------------------------
function check_epipoles_in_image( e1, e2, rows, cols )
% Check whether given epipoles are in the image or not
if abs( e1(3) ) < 1e-6 & abs( e2(3) ) < 1e-6, return; end;
e1 = e1 / e1(3);
e2 = e2 / e2(3);
if ( e1(1) <= cols & e1(1) >= 1 & e1(2) <= rows & e1(2) >= 1 ) | ...
( e2(1) <= cols & e2(1) >= 1 & e2(2) <= rows & e2(2) >= 1 )
err_msg = sprintf( 'epipole (%g,%g) or (%g,%g) is inside image', ...
e1(1:2), e2(1:2) );
error( [ err_msg, ' -- homography does not work in this case!' ] );
end;
%-----------------------------------------------------------------------------
function [P,Pprime] = get_canonical_cameras( F )
% Get the "canonical" cameras for given fundamental matrix
% according to Hartley and Zisserman (2004), p256, Result 9.14
% But ensure that the left 3x3 submatrix of Pprime is nonsingular
% using Result 9.15, that the general form is
% [ skewsym( e12 ) * F + e12 * v', k * e12 ] where v is an arbitrary
% 3-vector and k is an arbitrary scalar
P = [ 1 0 0 0
0 1 0 0
0 0 1 0 ];
e12 = null( F' );
M = skew( e12 ) * F + e12 * [1 1 1];
Pprime = [ M, e12 ];
%-----------------------------------------------------------------------------
function test_p_f( P, Pprime, F )
% Test that camera matrices Pprime and P are consistent with
% fundamental matrix F
% Meaning (Pprime*X)' * F * (P*X) = 0, for all X in 3space
% Get the epipole in camera 1 for camera 2
C2 = null( P );
eprime = Pprime * C2;
% Construct F from Pprime, P, and eprime
Fhat = skew( eprime ) * Pprime * pinv( P );
% Check that it's close to F
alpha = Fhat(:)\F(:);
if norm( alpha*Fhat-F ) > 1e-10
fprintf( 1, 'Warning: supplied camera matrices are inconsistent with F\n' );
else
fprintf( 1, 'Supplied camera matrices OK\n' );
end;
%-----------------------------------------------------------------------------
function H = map_to_infinity( x, cx, cy )
% Given a point and the desired origin (point of minimum projective
% distortion), compute a homograph H = G*R*T taking the point to the
% origin, rotating it to align with the X axis, then mapping it to
% infinity.
% First map cx,cy to the origin
T = [ 1 0 -cx
0 1 -cy
0 0 1 ];
x = T * x;
% Now rotate the translated x to align with the X axis.
cur_angle = atan2( x(2), x(1) );
R = [ cos( -cur_angle ), -sin( -cur_angle ), 0
sin( -cur_angle ), cos( -cur_angle ), 0
0, 0, 1 ];
x = R * x;
% Now the transformation G mapping x to infinity
if abs( x(3)/norm(x) ) < 1e-6
% It's already at infinity
G = eye(3)
else
f = x(1)/x(3);
G = [ 1 0 0
0 1 0
-1/f 0 1 ];
end;
H = G*R*T; |
idisp | image display tool
%
% IDISP(IM, OPTIONS) displays an image and allows interactive investigation
% such as pixel value, cross-section, histogram and zooming. The image is
% displayed in a figure with a toolbar across the top.
%
% IDISP(C, OPTIONS) horizontally concatenates the images in the cell array C
% and displays them as above.
%
% User interface:
%
% - Left clicking on a pixel will display its value in a box at the top.
%
% - The "line" button allows two points to be specified and a new figure
% displays intensity along a line between those points.
%
% - The "histo" button displays a histogram of the pixel values in a new figure.
% If the image is zoomed, the histogram is computed over only those pixels in
% view.
%
% The "zoom" button requires a left-click and drag to specify a box which
% defines the zoomed view.
%
% Options::
% 'ncolors',N number of colors in the color map (default 256)
% 'nogui' display the image without the GUI
% 'noaxes' no axes on the image
% 'noframe' no axes or frame on the image
% 'plain' no axes, frame or GUI
% 'bar' add a color bar to the image
% 'print',F write the image to file F in EPS format
% 'square' display aspect ratio so that pixels are squate
% 'wide' make figure very wide, useful for displaying stereo pair
% 'flatten' display image planes as 1G
% 'ynormal' y-axis increases upward, image is inverted
% 'cscale',C C is a 2-vector that specifies the grey value range that spans
% the colormap.
% 'xydata',XY XY is a cell array whose elements are vectors that span the
% x- and y-axes respectively.
% 'grey' color map: greyscale unsigned
% 'invert' color map: greyscale unsigned, zero is white, maximum value
% is black
% 'signed' color map: greyscale signed, positive is blue, negative is red,
% zero is black
% 'invsigned' color map: greyscale signed, positive is blue, negative is red,
% zero is white
% 'random' color map: random values, highlights fine structure
% 'dark' color map: greyscale unsigned, darker than 'grey', good for
% superimposed graphics
%
% Notes::
% - color images are displayed in true color mode: pixel triples map to display
% pixels
% - grey scale images are displayed in indexed mode: the image pixel value is
% mapped through the color map to determine the display pixel value.
% - the minimum and maximum image values are mapped to the first and last
% element of the color map, which by default ('greyscale') is the range black
% to white.
%
% See also IMAGE, CAXIS, COLORMAP, ICONCAT.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function idisp(im, varargin)
opt.ncolors = 256;
opt.gui = true;
opt.axes = true;
opt.square = false;
opt.wide = false;
opt.colormap = [];
opt.print = [];
opt.bar = false;
opt.frame = true;
opt.ynormal = false;
opt.cscale = [];
opt.xydata = [];
opt.plain = false;
opt.flatten = false;
opt.colormap_std = {[], 'grey', 'signed', 'invsigned', 'random', 'invert', 'dark'};
[opt,arglist] = tb_optparse(opt, varargin);
if opt.plain
opt.frame = false;
opt.gui = false;
end
if ~isempty(opt.print)
opt.gui = false;
end
if opt.flatten
im = reshape( im, size(im,1), size(im,2)*size(im,3) );
end
if length(arglist) ~= 0
warning(['Unknown options: ', arglist]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% command line invocation, display the image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display the image
clf
ud = [];
if iscell(im)
% image is a cell array
[im,ud.u0] = iconcat(im);
end
ud.size = size(im);
set(gca, 'CLimMode', 'Manual');
set(gca, 'CLim', [min(im(:)) max(im(:))]);
if ~isempty(opt.xydata)
hi = image(opt.xydata{1}, opt.xydata{2}, im);
else
hi = image(im, 'CDataMapping', 'scaled');
end
if opt.wide
set(gcf, 'units', 'norm');
pos = get(gcf, 'pos');
set(gcf, 'pos', [0.0 pos(2) 1.0 pos(4)]);
end
if isempty(opt.colormap_std)
if isempty(opt.colormap)
% default colormap
disp('default color map');
cmap = gray(opt.ncolors);
else
% load a Matlab color map
disp('matlab color map');
cmap = feval(opt.colormap);
end
else
% a builtin shorthand color map was specified
disp(['builtin color map: ', opt.colormap_std]);
switch opt.colormap_std
case 'random'
cmap = rand(opt.ncolors,3);
case 'dark'
cmap = gray(opt.ncolors)*0.5;
case 'grey'
cmap = gray(opt.ncolors);
case 'invert'
% invert the monochrome color map: black <-> white
cmap = gray(opt.ncolors);
cmap = cmap(end:-1:1,:);
case {'signed', 'invsigned'}
% signed color map, red is negative, blue is positive, zero is black
% inverse signed color map, red is negative, blue is positive, zero is white
cmap = zeros(opt.ncolors, 3);
opt.ncolors = bitor(opt.ncolors, 1); % ensure it's odd
ncm2 = ceil(opt.ncolors/2);
if strcmp(opt.colormap, 'signed')
% signed color map, red is negative, blue is positive, zero is black
for i=1:opt.ncolors
if i > ncm2
cmap(i,:) = [0 0 1] * (i-ncm2) / ncm2;
else
cmap(i,:) = [1 0 0] * (ncm2-i) / ncm2;
end
end
else
% inverse signed color map, red is negative, blue is positive, zero is white
for i=1:opt.ncolors
if i > ncm2
s = (i-ncm2)/ncm2;
cmap(i,:) = [1-s 1-s 1];
else
s = (ncm2-i)/ncm2;
cmap(i,:) = [1 1-s 1-s];
end
end
end
mn = min(im(:));
mx = max(im(:));
set(gca, 'CLimMode', 'Manual');
if mn < 0 && mx > 0
a = max(-mn, mx);
set(gca, 'CLim', [-a a]);
elseif mn > 0
set(gca, 'CLim', [-mx mx]);
elseif mx < 0
set(gca, 'CLim', [-mn mn]);
end
end
end
colormap(cmap);
if opt.bar
colorbar
end
if opt.frame
if opt.axes
xlabel('u (pixels)');
ylabel('v (pixels)');
else
set(gca, 'Xtick', [], 'Ytick', []);
end
else
set(gca, 'Visible', 'off');
end
if opt.square
set(gca, 'DataAspectRatio', [1 1 1]);
end
if opt.ynormal
set(gca, 'YDir', 'normal');
end
set(hi, 'CDataMapping', 'scaled');
if ~isempty(opt.cscale)
set(gca, 'Clim', opt.cscale);
end
figure(gcf); % bring to top
if opt.print
print(opt.print, '-depsc');
return
end
if opt.gui
set(gcf, 'MenuBar', 'none');
set(gcf, 'ToolBar', 'none');
htf = uicontrol(gcf, ...
'style', 'text', ...
'units', 'norm', ...
'pos', [.5 .935 .48 .05], ...
'background', [1 1 1], ...
'HorizontalAlignment', 'left', ...
'string', ' Machine Vision Toolbox for Matlab ' ...
);
ud.axis = gca;
ud.panel = htf;
ud.image = hi;
set(gca, 'UserData', ud);
set(hi, 'UserData', ud);
% show the variable name in the figure's title bar
varname = inputname(1);
if isempty(varname)
set(gcf, 'name', 'idisp');
else
set(gcf, 'name', sprintf('idisp: %s', varname));
end
% create pushbuttons
uicontrol(gcf,'Style','Push', ...
'String','line', ...
'Foregroundcolor', [0 0 1], ...
'Units','norm','pos',[0 .93 .1 .07], ...
'UserData', ud, ...
'Callback', @(src,event) idisp_callback('line', src) );
uicontrol(gcf,'Style','Push', ...
'String','histo', ...
'Foregroundcolor', [0 0 1], ...
'Units','norm','pos',[0.1 .93 .1 .07], ...
'UserData', ud, ...
'Callback', @(src,event) idisp_callback('histo', src) );
uicontrol(gcf,'Style','Push', ...
'String','zoom', ...
'Foregroundcolor', [0 0 1], ...
'Units','norm','pos',[.2 .93 .1 .07], ...
'Userdata', ud, ...
'Callback', @(src,event) idisp_callback('zoom', src) );
uicontrol(gcf,'Style','Push', ...
'String','unzoom', ...
'Foregroundcolor', [0 0 1], ...
'Units','norm','pos',[.3 .93 .15 .07], ...
'Userdata', ud, ...
'Callback', @(src,event) idisp_callback('unzoom', src) );
%'DeleteFcn', 'idisp(''cleanup'')', ...
set(gcf, 'Color', [0.8 0.8 0.9], ...
'WindowButtonDownFcn', @(src,event) idisp_callback('down', src), ...
'WindowButtonUpFcn', @(src,event) idisp_callback('up', src) );
% htf = uicontrol(gcf, ...
% 'style', 'text', ...
% 'units', 'norm', ...
% 'pos', [.6 0 .4 .04], ...
% 'ForegroundColor', [0 0 1], ...
% 'BackgroundColor', get(gcf, 'Color'), ...
% 'HorizontalAlignment', 'right', ...
% 'string', 'Machine Vision Toolbox for Matlab ' ...
% );
set(hi, 'UserData', ud);
end
set(hi, 'DeleteFcn', @(src,event) idisp_callback('idelete', src) );
set(gca, ...
'DeleteFcn', @(src,event) idisp_callback('destroy', src), ...
'NextPlot', 'replace', ...
'UserData', ud);
end
% invoked on a GUI event
function idisp_callback(cmd, src)
%disp(['in callback: ', cmd]);
if isempty(cmd)
% mouse push or motion request
h = get(gcf, 'CurrentObject'); % image
ud = get(h, 'UserData'); % axis
if ~isempty(ud)
cp = get(ud.axis, 'CurrentPoint');
x = round(cp(1,1));
y = round(cp(1,2));
try
imdata = get(ud.image, 'CData');
set(ud.panel, 'String', sprintf(' (%d, %d) = %s', x, y, num2str(imdata(y,x,:), 4)));
drawnow
end
end
else
switch cmd
case {'destroy', 'idelete'}
%fprintf('cleaning up figure\n');
clf
set(gcf, 'MenuBar', 'figure');
set(gcf, 'ToolBar', 'figure');
set(gcf, 'WindowButtonUpFcn', '');
set(gcf, 'WindowButtonDownFcn', '');
case 'cleanup'
%fprintf('cleaning up handlers\n');
set(gcf, 'WindowButtonDownFcn', '');
set(gcf, 'WindowButtonUpFcn', '');
case 'down',
% install pixel value inspector
set(gcf, 'WindowButtonMotionFcn', @(src,event) idisp_callback([], src) );
idisp_callback([], src);
case 'up',
set(gcf, 'WindowButtonMotionFcn', '');
case 'line',
h = get(gcf, 'CurrentObject'); % push button
ud = get(h, 'UserData');
set(ud.panel, 'String', 'Click first point');
[x1,y1] = ginput(1);
x1 = round(x1); y1 = round(y1);
set(ud.panel, 'String', 'Click last point');
[x2,y2] = ginput(1);
x2 = round(x2); y2 = round(y2);
set(ud.panel, 'String', '');
imdata = get(ud.image, 'CData');
dx = x2-x1; dy = y2-y1;
if abs(dx) > abs(dy),
x = x1:x2;
y = round(dy/dx * (x-x1) + y1);
figure
if size(imdata,3) > 1
set(gca, 'ColorOrder', eye(3,3), 'Nextplot', 'replacechildren');
n = size(imdata,1)*size(imdata,2);
z = [];
for i=1:size(imdata,3)
z = [z imdata(y+x*numrows(imdata)+(i-1)*n)'];
end
plot(x', z);
else
plot(imdata(y+x*numrows(imdata)))
end
else
y = y1:y2;
x = round(dx/dy * (y-y1) + x1);
figure
if size(imdata,3) > 1
set(gca, 'ColorOrder', eye(3,3), 'Nextplot', 'replacechildren');
n = size(imdata,1)*size(imdata,2);
z = [];
for i=1:size(imdata,3)
z = [z imdata(y+x*numrows(imdata)+(i-1)*n)'];
end
plot(z, y');
else
plot(imdata(y+x*numrows(imdata)))
end
end
title(sprintf('(%d,%d) to (%d,%d)', x1, y1, x2, y2));
xlabel('distance (pixels)')
ylabel('greyscale');
grid on
case 'histo',
h = get(gcf, 'CurrentObject'); % push button
ud = get(h, 'UserData');
imdata = get(ud.image, 'CData');
b = floor(axis); % bounds of displayed image
if b(1) == 0,
b = [1 b(2) 1 b(4)];
end
figure
imdata = double(imdata(b(3):b(4), b(1):b(2),:));
ihist(imdata);
case 'zoom',
% get the rubber band box
waitforbuttonpress
cp0 = floor( get(gca, 'CurrentPoint') );
rect = rbbox; % return on up click
cp1 = floor( get(gca, 'CurrentPoint') );
% determine the bounds of the ROI
top = cp0(1,2);
left = cp0(1,1);
bot = cp1(1,2);
right = cp1(1,1);
if bot |
idisplabel | Display an image with mask
%
% IDISPLABEL(IM, LABELIMAGE, LABELS) displays only those image pixels which
% belong to a specific class. IM is a greyscale NxM or color NxMx3 image, and
% LABELIMAGE is an NxM image containing integer pixel class labels for the
% corresponding pixels in IM. The pixel classes to be displayed are given by
% the elements of LABELS which is a scalar a vector of class labels.
% Non-selected pixels are displayed as white.
%
% IDISPLABEL(IM, LABELIMAGE, LABELS, BG) as above but the grey level of the
% non-selected pixels is specified by BG in the range 0 to 1.
%
% See also IBLOBS, ICOLORIZE, COLORSEG.
function idisplabel(im, label, select, bg)
if isscalar(select)
mask = label == select;
else
mask = zeros(size(label));
for s=select(:)',
mask = mask | (label == s);
end
end
if nargin < 4
bg = 1;
end
if ndims(im) == 3
mask = cat(3, mask, mask, mask);
end
im(~mask) = bg;
idisp(im, 'nogui');
shg |
iread | Read an image from file
%
% IM = IREAD() presents a file selection GUI from which the user can select
% an image file which is returned as 2D or 3D matrix. On subsequent calls
% the initial folder is as set on the last call.
%
% IM = IREAD(FILE, OPTIONS) reads the specified file and returns a matrix. If
% the path is relative it is searched for on Matlab search path.
% Wildcards are allowed in file names. If multiple files match a 3D or 4D image
% is returned where the last dimension is the number of images in the sequence.
%
% Options::
% 'uint8' return an image with 8-bit unsigned integer pixels in
% the range 0 to 255
% 'single' return an image with single precision floating point pixels
% in the range 0 to 1.
% 'double' return an image with double precision floating point pixels
% in the range 0 to 1.
% 'grey' convert image to greyscale if it's color using ITU rec 601
% 'grey_601' ditto
% 'grey_709' convert image to greyscale if it's color using ITU rec 709
% 'gamma',G gamma value, either numeric or 'sRGB'
% 'reduce',R decimate image by R in both dimensions
% 'roi',R apply the region of interest R to each image,
% where R=[umin umax; vmin vmax].
% Notes::
% - a greyscale image is returned as an NxM matrix
% - a color image is returned as an NxMx3 matrix
% - a greyscale image sequence is returned as an NxMxP matrix where P is the sequence length
% - a color image sequence is returned as an NxMx3xP matrix where P is the sequence length
%
% See also IDISP, IMONO, IGAMMA, IMWRITE, PATH.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function [I,info] = iread(filename, varargin)
persistent mypath
% options
%
% 'float
% 'uint8
% 'grey'
% 'gray'
% 'grey_601'
% 'grey_709'
% 'grey_value'
% 'gray_601'
% 'gray_709'
% 'gray_value'
% 'reduce', n
opt.type = {[], 'double', 'single', 'uint8'};
opt.mkGrey = {[], 'grey', 'gray', 'mono', '601', 'grey_601', 'grey_709'};
opt.gamma = [];
opt.reduce = 1;
opt.roi = [];
opt.disp = [];
opt = tb_optparse(opt, varargin);
im = [];
if nargin == 0,
% invoke file browser GUI
[file, npath] = uigetfile(...
{'*.png;*.pgm;*.ppm;*.jpg;*.tif', 'All images';
'*.pgm', 'PGM images';
'*.jpg', 'JPEG images';
'*.gif;*.png;*.jpg', 'web images';
'*.*', 'All files';
}, 'iread');
if file == 0,
fprintf('iread canceled from GUI\n');
return; % cancel button pushed
else
% save the path away for next time
mypath = npath;
filename = fullfile(mypath, file);
im = loadimg(filename, opt);
end
elseif (nargin == 1) & exist(filename,'dir'),
% invoke file browser GUI
if isempty(findstr(filename, '*')),
filename = strcat(filename, '/*.*');
end
[file,npath] = uigetfile(filename, 'iread');
if file == 0,
fprintf('iread canceled from GUI\n');
return; % cancel button pushed
else
% save the path away for next time
mypath = npath;
filename = fullfile(mypath, file);
im = loadimg(filename, opt);
end
else
% some kind of filespec has been given
if ~isempty(strfind(filename, '*')) | ~isempty(strfind(filename, '?')),
% wild card files, eg. 'seq/*.png', we need to look for a folder
% seq somewhere along the path.
[pth,name,ext] = fileparts(filename);
if opt.verbose
fprintf('wildcard lookup: %s %s %s\n', pth, name, ext);
end
% search for the folder name along the path
folderonpath = pth;
for p=path2cell(path)' % was userpath
if exist( fullfile(p{1}, pth) ) == 7
folderonpath = fullfile(p{1}, pth);
if opt.verbose
fprintf('folder found\n');
end
break;
end
end
s = dir( fullfile(folderonpath, [name, ext])); % do a wildcard lookup
if length(s) == 0,
error('no matching files found');
end
for i=1:length(s)
im1 = loadimg( fullfile(folderonpath, s(i).name), opt);
if i==1
% preallocate storage, much quicker
im = zeros([size(im1) length(s)], class(im1));
end
if ndims(im1) == 2
im(:,:,i) = im1;
elseif ndims(im1) == 3
im(:,:,:,i) = im1;
end
end
else
% simple file, no wildcard
if strncmp(filename, 'http://', 7)
im = loadimg(filename, opt);
elseif exist(filename)
im = loadimg(filename, opt);
else
% see if it exists on the Matlab search path
for p=path2cell(userpath)'
if exist( fullfile(p{1}, filename) ) > 0
im = loadimg(fullfile(p{1}, filename), opt);
break;
end
end
end
end
end
if isempty(im)
error(sprintf('cant open file: %s', filename));
end
if nargout > 0
I = im;
if nargout > 1
info = imfinfo(filename);
end
else
% if no output arguments display the image
if ndims(I) <= 3
idisp(I);
end
end
end
function im = loadimg(name, opt)
% now we read the image
im = imread(name);
if opt.verbose
if ndims(im) == 2
fprintf('loaded %s, %dx%d\n', name, size(im,2), size(im,1));
elseif ndims(im) == 3
fprintf('loaded %s, %dx%dx%d\n', name, size(im,2), size(im,1), size(im,3));
end
end
% optionally convert it to greyscale using specified method
if ~isempty(opt.mkGrey) && (ndims(im) == 3)
im = imono(im, opt.mkGrey);
end
% optionally chop out a roi
if ~isempty(opt.roi)
im = iroi(im, opt.roi);
end
% optionally decimate it
if opt.reduce > 1,
im = im(1:opt.reduce:end, 1:opt.reduce:end, :);
end
% optionally convert to specified numeric type
if ~isempty(opt.type)
if isempty(findstr(opt.type, 'int'))
im = cast(im, opt.type) / double(intmax(class(im)));
else
im = cast(im, opt.type);
end
end
% optionally gamma correct it
if ~isempty(opt.gamma)
im = igamma(im, opt.gamma);
end
if opt.disp
idisp(im);
end
end |
pgmfilt | Pipe image through PGM utility
%
% OUT = PGMFILT(IM, PGMCMD) pipes the image IM through an Unix filter program
% and returns its output as an image. The program given by the string PGMCMD
% must accept and return images in PGM format.
%
% See also PNMFILT, IREAD.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function im2 = pgmfilt(im, cmd)
% MATLAB doesn't support pipes, so it all has to be done via temp files...
% make up two file names
fname = tempname;
fname2 = tempname;
imwrite(im, fname, 'pgm');
%cmd
unix([cmd ' < ' fname ' > ' fname2]);
%fname2
im2 = imread(fname2, 'pgm');
unix(['/bin/rm ' fname ' ' fname2]); |
pnmfilt | Pipe image through PNM utility
%
% OUT = PNMFILT(IM, PNMCMD) pipes the image IM through an Unix filter program
% and returns its output as an image. The program given by the string PNMCMD
% must accept and return images in PNM format.
%
% See also PGMFILT, IREAD.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function im2 = pnmfilt(im, cmd)
% MATLAB doesn't support pipes, so it all has to be done via
% temp files :-(
% make up two file names
ifile = sprintf('%s.pnm', tempname);
ofile = sprintf('%s.pnm', tempname);
imwrite(im, ifile, 'pgm');
%cmd
unix([cmd ' < ' ifile ' > ' ofile]);
im2 = double( imread(ofile) );
unix(['/bin/rm ' ifile ' ' ofile]); |
bresenham | Create points on line
%
% P = BRESENHAM(X1, Y1, X2, Y2) is a list of integer coordinates for
% points lying on the line segement (X1,Y1) to (X2,Y2). Endpoints
% must be integer.
%
% P = BRESENHAM(P1, P2) as above but P1=[X1,Y1] and P2=[X2,Y2].
%
% See also ICANVAS.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function p = bresenham(x1, y1, x2, y2)
if nargin == 2
p1 = x1; p2 = y1;
x1 = p1(1); y1 = p1(2);
x2 = p2(1); y2 = p2(2);
elseif nargin ~= 4
error('expecting 2 or 4 arguments');
end
x = x1;
if x2 > x1
xd = x2-x1;
dx = 1;
else
xd = x1-x2;
dx = -1;
end
y = y1;
if y2 > y1
yd = y2-y1;
dy = 1;
else
yd = y1-y2;
dy = -1;
end
p = [];
if xd > yd
a = 2*yd;
b = a - xd;
c = b - xd;
while 1
p = [p; x y];
if all([x-x2 y-y2] == 0)
break
end
if b < 0
b = b+a;
x = x+dx;
else
b = b+c;
x = x+dx; y = y+dy;
end
end
else
a = 2*xd;
b = a - yd;
c = b - yd;
while 1
p = [p; x y];
if all([x-x2 y-y2] == 0)
break
end
if b < 0
b = b+a;
y = y+dy;
else
b = b+c;
x = x+dx; y = y+dy;
end
end
end
end |
closest | Find matching points in N-dimensional space.
%
% K = CLOSEST(A, B) returns point correspondence for N-dimensional points. A
% is NxNA and B is NxNB. K is 1 x NA and the element J = K(I) indicates that
% the I'th column of A is closest to the Jth column of B. That is, A(:,I)
% is closest to B(:,J).
%
% Notes::
% - is a MEX file. |
col2im | Convert pixel vector to image
%
% OUT = COL2IM(PIX, IMSIZE) is an image comprising the pixel values in PIX which
% is an LxP matrix with one row per pixel. The result is an NxMxP image
% where L = NxM. IMSIZE is a 2-vector (N,M).
%
% OUT = COL2IM(PIX, IM) as above but the dimensions of OUT are the same as IM.
%
% See also IM2COL.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function im = col2im(col, img)
ncols = numcols(col);
if ~((ncols == 1) | (ncols == 3)),
error('must be 1 or 3 columns');
end
if length(img) == 2,
sz = img;
elseif length(img) == 3,
sz = img(1:2);
else
sz = size(img);
sz = sz(1:2);
end
if ncols == 3,
sz = [sz 3];
end
sz
im = reshape(col, sz); |
distance | DISTANCE Euclidean distances between sets of points
%
% D = DISTANCE(A,B) is the Euclidean distances between L-dimensional points
% described by the matrices A (LxM) and B (LxN) respectively. The distance
% D is MxN and element D(I,J) is the distance between points A(I) and D(J).
%
% Notes::
% - This fully vectorized (VERY FAST!)
% - It computes the Euclidean distance between two vectors by:
%
% ||A-B|| = sqrt ( ||A||^2 + ||B||^2 - 2*A.B )
%
% Example::
% A = rand(400,100); B = rand(400,200);
% d = distance(A,B);
%
% Author::
% Roland Bunschoten,
% University of Amsterdam,
% Intelligent Autonomous Systems (IAS) group,
% Kruislaan 403 1098 SJ Amsterdam,
% tel.(+31)20-5257524,
% bunschot@wins.uva.nl
% Last Rev : Oct 29 16:35:48 MET DST 1999
% Tested : PC Matlab v5.2 and Solaris Matlab v5.3
% Thanx : Nikos Vlassis
% Copyright notice: You are free to modify, extend and distribute
% this code granted that the author of the original code is
% mentioned as the original author of the code.
function d = distance(a,b)
if (nargin ~= 2)
error('Not enough input arguments');
end
if (size(a,1) ~= size(b,1))
error('A and B should be of same dimensionality');
end
aa=sum(a.*a,1); bb=sum(b.*b,1); ab=a'*b;
d = sqrt(abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab)); |
im2col | Convert an image to pixel per row format
%
% OUT = IM2COL(IM) returns the image as a pixel vector. Each row is a pixel
% value which is a P-vector where the image in NxMxP. The pixels are in
% image column order and there are NxM rows.
%
% See also COL2IM.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function c = im2col(im)
im = shiftdim(im, 2);
if ndims(im) == 3,
c = reshape(im, 3, [])';
else
c = reshape(im, 1, [])';
end |
imeshgrid | Domain matrices for image
%
% [U,V] = IMESHGRID(IM) return matrices that describe the domain of image IM
% and can be used for the evaluation of functions over the image. The
% element U(v,u) = u and V(v,u) = v.
%
% [U,V] = IMESHGRID(W, H) as above but the domain is WxH.
%
% [U,V] = IMESHGRID(SIZE) as above but the domain is described size which is
% scalar SIZExSIZE or a 2-vector [W H].
%
% See also MESHGRID.
function [U,V] = imeshgrid(a1, a2)
if nargin == 1
if length(a1) == 1
% we specified a size for a square output image
[U,V] = meshgrid(1:a1, 1:a1);
elseif length(a1) == 2
% we specified a size for a rectangular output image (w,h)
[U,V] = meshgrid(1:a1(1), 1:a1(2));
elseif ndims(a1) >= 2
[U,V] = meshgrid(1:numcols(a1), 1:numrows(a1));
else
error('incorrect argument');
end
elseif nargin == 2
[U,V] = meshgrid(1:a1, 1:a2);
end |
iscolor | Test if argument is a color image
%
% ISCOLOR(IM) is true (1) if IM is a color image, that is, it has 3 or
% more dimensions.
function s = iscolor(im)
s = isnumeric(im) && size(im,1) > 1 && size(im,2) > 1 && size(im,3) == 3; |
isize | Return size of image
%
% N = ISIZE(IM,D) returns the size of the D'th dimension of IM.
%
% [H,W] = ISIZE(IM) returns the image height H and width W.
%
% [H,W,P] = ISIZE(IM) returns the image height H and width W and number of
% planes.
%
% See also SIZE.
function [o1,o2,o3] = isize(im, idx)
if nargin == 2
o1 = size(im, idx);
else
s = size(im);
o1 = s(2); % width, number of columns
if nargout > 1
o2 = s(1); % height, number of rows
end
if nargout > 2
if ndims(im) == 2
o3 = 1;
else
o3 = s(3);
end
end
end |
kmeans | K-means clustering
%
% [L,C] = KMEANS(X, K, OPTIONS) performs K-means clustering for multi-dimensional data
% points X (NxD) where N is the number of points, and D is the dimension. K is
% the number of clusters. L is a vector (Nx1) whose elements indicates which
% cluster the corresponding element of X belongs to. The optional return
% value C (KxD) contains the cluster centroids.
%
% [L,C] = KMEANS(X, K, C0) as above but the initial clusters C0 (KxD) is given
% and row I is the initial estimate of the centre of cluster I.
%
% L = KMEANS(X, C) similar to above but the clustering step is not performed,
% it is assumed to have been completed previously. C (KxD) contains the cluster
% centroids and L (Nx1) indicates which cluster the corresponding element of X
% is closest to.
%
% Options::
% 'random' randomly choose K points from X
% 'spread' randomly choose K values within the hypercube spanned
% by X.
%
% Reference::
% Tou and Gonzalez,
% Pattern Recognition Principles,
% pp 94
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
% add Kaufman initialization
%
% 1) Initialize c1 with the most centrally located input sample.
% 2) For i = 2, . . . , K, each ci is initialized in the following
% way: for each non-selected input sample xj , calculate its
% summed distance to the other non-selected input samples xl ,
% who are closer to xj than to their respective nearest seed
% clusters,
% http://www.comp.nus.edu.sg/~tancl/Papers/IJCNN04/he04ijcnn.pdf
% http://citeseer.ist.psu.edu/398385.html
function [label,centroid,resid] = kmeans(x, K, varargin)
% TODO update doco for assignment mode
% option to loop N times and take lowest residual
% return residual from assignment mode
deb = 0;
n = numcols(x);
if numcols(K) > 1 && numrows(x) == numrows(K)
% kmeans(x, centres)
% then return closest clusters
label = closest(x, K);
if nargout > 1
centroid = K;
end
if nargout > 2
resid = norm( x - K(:,label) );
end
return
end
opt.plot = false;
opt.init = {'random', 'spread'};
[opt,z0] = tb_optparse(opt, varargin);
if opt.plot && numrow(x) > 3
warning('cant plot for more than 3D data');
opt.plot = false;
end
if ~isempty(z0)
% an initial condition was supplied
if numcols(z0) ~= K,
error('initial cluster length should be k');
end
if numrows(z0) ~= numrows(x),
error('number of dimensions of z0 must match dimensions of x');
end
else
% determine initial condition
if strcmp(opt.init, 'random')
% select K points from the set given as initial cluster centres
k = randi(n, K, 1);
z0 = x(:,k);
elseif strcmp(opt.init, 'spread')
% select K points from within the hypercube defined by the points
mx = max(x')';
mn = min(x')';
z0 = (mx-mn) * rand(1,K) + mn * ones(1,K);
end
end
% z is the centroid
% zp is the previous centroid
% s is the vector of cluster labels corresponding to rows in x
z = z0;
%
% step 1
%
zp = z; % previous centroids
s = zeros(1,n);
sp = s;
iterating = 1;
k = 1;
iter = 0;
while iterating,
iter = iter + 1;
tic
t0 = toc;
%
% step 2
%
s = closest(x, z);
%
% step 3
%
for j=1:K
k = find(s==j);
if isempty(k)
if strcmp(opt.init, 'random')
k = randi(n, 1, 1);
z0 = x(:,k);
elseif strcmp(opt.init, 'spread')
% this cluster has no elements, randomly assign it
zp(:,j) = (mx-mn) * rand(1,1) + mn;
end
else
% else, assign it to the mean of its elements
zp(:,j) = mean( x(:,k)' );
dd = sum( (repmat(zp(:,j),1,length(k)) - x(:,k)).^2 );
dd = dd / numrows(x);
maxerr(j) = max( sqrt(dd) );
end
end
%
% step 4
%
% determine the change in cluster centres over the last step
delta = sum( sum( (z - zp).^2 ) );
if opt.verbose
t = toc;
fprintf('%d: norm=%g, delta=%g, took %.1f seconds\n', iter, mean(maxerr), delta, t);
end
if delta == 0,
fprintf('delta to zero\n');
iterating = 0;
end
z = zp;
if opt.plot
if numrows(z) == 2
plot(z(1,:), z(2,:));
else
plot3(z(1,:), z(2,:), z(3,:));
end
pause(.1);
end
% if no point assignments changed then we are done
if all(sp == s)
fprintf('no point assignments changed\n');
iterating = 0;
end
sp = s;
end
if nargout == 0
% if no output arguments display results
if numrows(z) < 5
for i=1:K,
fprintf('cluster %d: %s (%d elements)\n', i, ...
sprintf('%11.4g ', z(i,:)), length(find(s==i)));
end
end
end
if nargout > 0
centroid = z;
end
if nargout > 1
label = s;
end
if nargout > 2
% compute the residual
resid = norm( x - z(:,s) );
end
dt = toc - t0;
if dt > 10
fprintf('that took %.1f seconds\n', dt);
elseif dt > 60
fprintf('that took %.1f minutes\n', dt/60);
end |
medfilt1 | One-dimensional median filter
%
% Y = MEDFILT1(X, W) is the median filter of the signal X computed over a
% sliding window of width W.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function m = medfilt1(s, w)
if nargin == 1,
w = 5;
end
s = s(:)';
w2 = floor(w/2);
w = 2*w2 + 1;
n = length(s);
m = zeros(w,n+w-1);
s0 = s(1); sl = s(n);
for i=0:(w-1),
m(i+1,:) = [s0*ones(1,i) s sl*ones(1,w-i-1)];
end
m = median(m);
m = m(w2+1:w2+n); |
ransac | Random sample and consensus
%
% M = RANSAC(FUNC, X, T, OPTIONS) robustly fit the data X to the function FUNC
% using the RANSAC algorithm. X contains corresponding point data, one column
% per point pair. RANSAC determines the subset of points (inliers) that best
% fit the model described by the function FUNC. T is a threshold applied
% to the model FUNC applied to a data point, if too high the point is
% considered an outlier.
%
% [M,IN] = RANSAC(FUNC, X, T, OPTIONS) as above but returns the vector IN of
% column indices of X that describe the inlier point set.
%
% [M,IN,RESID] = RANSAC(FUNC, X, T, OPTIONS) as above but returns the final
% residual of applying FUNC to the inlier set.
%
% OUT = FUNC(R) is the function passed to RANSAC and it must accept
% a single argument R which is a structure:
%
% R.cmd string the operation to perform which is either
% R.debug logical display what's going on
% R.X 6xN data to work on, N point pairs
% R.t 1x1 threshold
% R.theta 3x3 estimated quantity to test
% R.misc cell private data
%
% The function return value is also a structure:
%
% OUT.s 1x1 sample size
% OUT.X 2DxN conditioned data
% OUT.misc cell private data
% OUT.inlier 1xM list of inliers
% OUT.valid logical if data is valid for estimation
% OUT.theta 3x3 estimated quantity
% OUT.resid 1x1 model fit residual
%
% The values of R.cmd are:
% 'size' return the minimum number of points required to compute
% an estimate to OUT.s
% 'condition' condition the point data to X -> OUT.X
% 'decondition' decondition the estimated model data to R.theta -> OUT.theta
% 'valid' true if a set of points is not degenerate, that is they
% will produce a model. This is used to discard random samples
% that do not result in useful models -> OUT.valid
% 'estimate' [OUT.theta,OUT.resid] = EST(R.X) returns the best fit model
% and residual for the subset of points R.X. If this function
% cannot fit a model then OUT.theta = []. If multiple models
% are found OUT.theta is a cell array.
% 'error' [OUT.inlier,OUT.theta] = ERR(R.theta,R.X,T) evaluates the
% distance from the model(s) R.theta to the points R.X and
% returns the best model OUT.theta and the subset of R.X that
% best supports (most inliers) that model.
%
% For some algorithms (eg. fundamental matrix) it is necessary to condition
% the data to improve the accuracy of model estimation. The data transform
% parameters are kept in the .misc element and used in the final step to
% transform the estimate based on conditioned data to original data, the
% deconditioning operation.
%
% Options::
% 'maxTrials',N maximum number of iterations (default 2000)
% 'maxDataTrials',N maximum number of attempts to select a non-degenerate
% data set (default 100)
%
% References::
% - M.A. Fishler and R.C. Boles. "Random sample concensus: A paradigm
% for model fitting with applications to image analysis and automated
% cartography". Comm. Assoc. Comp, Mach., Vol 24, No 6, pp 381-395, 1981
%
% - Richard Hartley and Andrew Zisserman. "Multiple View Geometry in
% Computer Vision". pp 101-113. Cambridge University Press, 2001
%
% Author::
% Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk at csse uwa edu au
% http://www.csse.uwa.edu.au/~pk
%
% See also FMATRIX, HOMOGRAPHY.
% Copyright (c) 2003-2006 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk at csse uwa edu au
% http://www.csse.uwa.edu.au/~pk
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
%
% May 2003 - Original version
% February 2004 - Tidied up.
% August 2005 - Specification of distfn changed to allow model fitter to
% return multiple models from which the best must be selected
% Sept 2006 - Random selection of data points changed to ensure duplicate
% points are not selected.
% February 2007 - Jordi Ferrer: Arranged warning printout.
% Allow maximum trials as optional parameters.
% Patch the problem when non-generated data
% set is not given in the first iteration.
% August 2008 - 'feedback' parameter restored to argument list and other
% breaks in code introduced in last update fixed.
% December 2008 - Octave compatibility mods
function [M, inliers, resid] = ransac(fittingfn, x, t, varargin);
%useRandomsample = ~exist('randsample');
opt.maxTrials = 2000;
opt.maxDataTrials = 100;
opt = tb_optparse(opt, varargin);
[rows, npts] = size(x);
p = 0.99; % Desired probability of choosing at least one sample
% free from outliers
bestM = NaN; % Sentinel value allowing detection of solution failure.
trialcount = 0;
bestscore = 0;
N = 1; % Dummy initialisation for number of trials.
out = invoke('size', fittingfn, opt.verbose, []);
s = out.s;
in.X = x;
out = invoke('condition', fittingfn, opt.verbose, in);
x = out.X;
misc = out.misc;
while N > trialcount
% Select at random s datapoints to form a trial model, M.
% In selecting these points we have to check that they are not in
% a degenerate configuration.
degenerate = 1;
count = 1;
while degenerate
% Generate s random indicies in the range 1..npts
% (If you do not have the statistics toolbox, or are using Octave,
% use the function RANDOMSAMPLE from my webpage)
if 0
if useRandomsample
ind = randomsample(npts, s);
else
ind = randsample(npts, s);
end
else
ind = randi(npts, 1, s);
end
% Test that these points are not a degenerate configuration.
in.X = x(:,ind);
out = invoke('valid', fittingfn, opt.verbose, in);
degenerate = ~out.valid;
if ~degenerate
% Fit model to this random selection of data points.
% Note that M may represent a set of models that fit the data in
% this case M will be a cell array of models
in.X = x(:,ind);
out = invoke('estimate', fittingfn, opt.verbose, in);
M = out.theta;
% Depending on your problem it might be that the only way you
% can determine whether a data set is degenerate or not is to
% try to fit a model and see if it succeeds. If it fails we
% reset degenerate to true.
if isempty(M)
degenerate = 1;
end
end
% Safeguard against being stuck in this loop forever
count = count + 1;
if count > opt.maxDataTrials
warning('Unable to select a nondegenerate data set');
break
end
end
% Once we are out here we should have some kind of model...
% Evaluate distances between points and model returning the indices
% of elements in x that are inliers. Additionally, if M is a cell
% array of possible models 'distfn' will return the model that has
% the most inliers. After this call M will be a non-cell object
% representing only one model.
in.t = t;
in.theta = M;
in.X = x;
out = invoke('error', fittingfn, opt.verbose, in);
inliers = out.inliers;
M = out.theta;
% Find the number of inliers to this model.
ninliers = length(inliers);
if ninliers > bestscore % Largest set of inliers so far...
bestscore = ninliers; % Record data for this model
bestinliers = inliers;
bestM = M;
% Update estimate of N, the number of trials to ensure we pick,
% with probability p, a data set with no outliers.
fracinliers = ninliers/npts;
pNoOutliers = 1 - fracinliers^s;
pNoOutliers = max(eps, pNoOutliers); % Avoid division by -Inf
pNoOutliers = min(1-eps, pNoOutliers);% Avoid division by 0.
N = log(1-p)/log(pNoOutliers);
end
trialcount = trialcount+1;
if opt.verbose > 1
fprintf('trial %d out of %d \r',trialcount, ceil(N));
end
% Safeguard against being stuck in this loop forever
if trialcount > opt.maxTrials
warning( ...
sprintf('ransac reached the maximum number of %d trials',...
opt.maxTrials));
break
end
end
fprintf('\n');
if ~isnan(bestM) % We got a solution
M = bestM;
inliers = bestinliers;
else
M = [];
inliers = [];
error('ransac was unable to find a useful solution');
end
% Now do a final least squares fit on the data points considered to
% be inliers.
[M,resid] = feval(fittingfn, x(:,inliers)); % with conditioned data
% then decondition it
in.theta = M;
in.misc = misc;
out = invoke('decondition', fittingfn, opt.verbose, in);
M = out.theta;
if opt.verbose || (nargout == 0)
fprintf('%d trials\n', trialcount);
fprintf('%d outliers\n', npts-length(inliers));
fprintf('%g final residual\n', resid);
end
end
function out = invoke(cmd, func, verbose, in)
if isempty(in) || (nargin < 4)
in = [];
end
in.debug = verbose > 1;
in.cmd = cmd;
if verbose > 2,
fprintf('---------------------------------------------\n');
in
end
out = feval(func, in);
if verbose > 2,
out
end
end |
zcross | Zero-crossing detector
%
% IZ = ZCROSS(IM) is a binary image with pixels set corresponding to
% zero crossings in the signed image IM. That is, pixels are set at positive
% pixels adajacent to a negative pixel.
%
% Notes::
% - can be used in association with a Lapalacian of Gaussian image to
% determine edges.
%
% See also ILOG.
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB. If not, see .
function iz = zcross(im)
z = zeros([size(im), 4]);
K = [1 1 0 ;1 1 0; 0 0 0];
z(:,:,1) = conv2(im, K, 'same');
K = [0 1 1 ;0 1 1; 0 0 0];
z(:,:,2) = conv2(im, K, 'same');
K = [0 0 0; 1 1 0 ;1 1 0];
z(:,:,3) = conv2(im, K, 'same');
K = [0 0 0; 0 1 1 ;0 1 1];
z(:,:,4) = conv2(im, K, 'same');
maxval = max(z, [], 3);
minval = min(z, [], 3);
iz = (maxval > 0) & (minval < 0); |