How to Extract High density pixels in a RGB Image?

How to Extract High density pixels in a RGB Image?

[Walter Roberson]

Please explain your question, with an example.

[harjeet singh]

do you mean high density by histogram?

[Image Analyst]

Not sure what you mean. In a radiograph (x-ray), the material density is indicated by the gray level. Do you have radiographs?

[Selva Karna]

My Need is how to skip or reduce glare light in a RGB image?

[Guillaume]

Not sure what a 'high density pixel' is or what that has to do with your latest comment.

Anyway, the best way to reduce glare light in an image is to actually change your imaging setup so that the glare light is not there in the first place.

Cleaning up mess is always more difficult than not generating the mess in the first place.

[Image Analyst]

Right, I remember this image and I suggested the same thing as Guillaume. Twice now since you just reposted it. So I'll make the suggestion a third time here. Please take our advice.

Answer 1

[John BG]

perhaps you mean high saturation pixels.

highly saturated colour is that with low or no contents of white.

White is flat colour spectrum. So the flatter is the RGB(or equivalent Cyan Magenta Yellow) base.

Red is [1 0 0], or in JPEG: [255 0 0] this is 100% saturated Red.

[1 .5 .5] or [255 127 127] would be a half saturated Red

[255 255 0] is a 100% saturated Yellow. [255 255 127] is 'half clean' Yellow.

Cyan [0 1 1] Mag [1 0 1]

Another possibility is that you are just referring to the amount of Luminance (in analog TV is known as Yr) that is just the amount common to all 3 base components.

Luminance is grey scaling, which is also the output of Infrared sensors, all the colouring of IR images comes afterwards to help 'read' what is hot and what is not.

Hope it helps

John [e-mail]

[Selva Karna]

Can you help me JOHAN sir ? i need how to count number object {GLASS } Above Images? Please refer me sir?

Answer 2

[John BG]

Hi Selva

1.- acquire image

A=imread('count_ribbons.jpg')
[im_1 im_2 im_3]=size(A) % [y x 3] or [vertical horizontal 3] or [rows columns 3]

to Image Analyst: let's work on the basis that there is no better image. We all see it's just a bunch of foil ribbons attached on a green cardboard, but think of it as a Radar or Sonar sample. It's what it is, right now you don't have anything better.

2.- check variance of RGB

surf(var(double(A),0,3))
surf(std(double(A),0,3))

3.- From var and std, there are heavily attenuated areas, especially on the upper half of the image.

Reading the image bottom up is as if target departing because the lines with ribbons get broader.

To process a bunch of lines without taking into account Doppler, let's take a strip made of the lowest 40 lines, process them, and answer the amount that clearly shows up as the mode among these 40 lines.

Before processing the strip, let's find out what processing is required for each line

4.- acquire a line within the last 40 lines

figure(1);imshow(A)
p3=ginput(2); p3=uint64(floor(p3)) % input 2 points excluding line % for instance:
% p3 = 33 394
% 403 38

5.- show the sample line

py_ref=min(p3(2,2),p3(1,2))
im_line3_1=A(py_ref,[min(p3(1,1),p3(2,1)):max(p3(1,1),p3(2,1))],1) % Red
im_line3_2=A(py_ref,[min(p3(1,1),p3(2,1)):max(p3(1,1),p3(2,1))],2) % Green
im_line3_3=A(py_ref,[min(p3(1,1),p3(2,1)):max(p3(1,1),p3(2,1))],3) % Blue
v3_ref=[1:1:(abs(double(p3(1,1))‐double(p3(2,1))))+1]
figure(2);plot(v3_ref,im_line3_1,'r',v3_ref,im_line3_2,'g',v3_ref,im_line3_3,'b')

observing the sample line, the ribbons are the valleys. The peaks are the green/white thin

6.- To count peaks coinciding with ribbons let's invert sample line:

im_line3_1=255‐im_line3_1
im_line3_2=255‐im_line3_2
im_line3_3=255‐im_line3_3

% figure(3);plot(v3_ref,im_line3_1,'r',v3_ref,im_line3_2,'g',v3_ref,im_line3_3,'b')

7.- Red component does not help that much, so carry on with Green and Blue only, (G+B) and amplify a bit

L1=3*(double(im_line3_2)+double(im_line3_3)) % amplify a bit
nL1=[1:1:length(L1)]
% figure(4);plot(L1)

8.- remove DC

L1=L1-mean(L1)
hold on;figure(4);plot(L1)

9.- 1st assault counting peaks, try different values and see how the amount of ribbons detected varies:

numel(findpeaks(L1,'minpeakheight',80))
numel(findpeaks(L1,'minpeakheight',100))
numel(findpeaks(L1,'minpeakheight',140))
numel(findpeaks(L1,'minpeakheight',170))
numel(findpeaks(L1,'minpeakheight',190))
sweep findpeaks parameter 'minpeakheight'
peaks_prob_1=0
for k=80:1:250
peaks_prob_1=[peaks_prob_1 numel(findpeaks(L1,'minpeakheight',k))];
end
figure(5);plot(peaks_prob_1);grid on
v_prob_1=[1:length(peaks_prob_1)]
figure(6)
hist(peaks_prob_1,v_prob_1) % misleading, the true amount of ribbons does not show

peaks count

peaks count histogram, cannot choose

10.- It's a really short sample, spectrum analysis needs as many samples as possible, and the data of interest to repeat, yet let's see DTFT:

    nL1=[1:length(L1)] % just copied lines from another script with DTFT
    x=L1;n=[1:1:length(L1)]
    x=double(x)
    k=0:500;w=(pi/500)*k
    X=x*(exp(‐j*pi/500)).^(n'*k)
    X=x*(exp(‐j*pi/500)).^(n'*k)
    magX=abs(X);angX=angle(X);realX=real(X);imagX=imag(X)
    figure(15);subplot(2,2,1);plot(k/500,magX);grid
    xlabel('freq[rad]');title('Magnitude(X)')
    figure(15);subplot(2,2,3);plot(k/500,angX/pi);grid
    xlabel('freq[rad]');title('Phase(X)')
    figure(15);subplot(2,2,2);plot(k/500,realX);grid
    xlabel('freq[rad]');title('Real(X)')
    figure(15);subplot(2,2,4);plot(k/500,imagX);grid
    xlabel('freq[pi]');title('Imag(X)')

count cycles:

371/2*.164=30.42 % do we round up? or down? error by 1 or 2 ribbons too many?

11.- and how well does the famous FFT?

    NFFT=2^nextpow2(length(L1)) % 512 freq points [0 2pi]
    absL1=abs(fft(L1,NFFT))
    figure(5);plot(absL1)

insert06

again count cycles:

371/2*43/256 % = 31.15 almost there but can't tell exactly

12.- try command FINDPEAKS and parameter 'minpeakprominence'

    [pks,locs]=findpeaks(L1,'minpeakprominence',10)
    [pks,locs]=findpeaks(L1,'minpeakprominence',15)
    [pks,locs]=findpeaks(L1,'minpeakprominence',20)
    [pks,locs]=findpeaks(L1,'minpeakprominence',25)
    [pks,locs]=findpeaks(L1,'minpeakprominence',50)
    [pks,locs]=findpeaks(L1,'minpeakprominence',60)
    % minpeakprominence=60 still one ribbon too many
    [pks,locs]=findpeaks(L1,'minpeakprominence',60);figure(7);plot(nL1,L1,locs,pks+4,'kv');grid numel(pks) %=31 and depending upon what line picked up

It is clear that more than one line has to be taken into account to make a decision.

13.- frequency contents with command PLOMB:

Fs=1;Fmax=1/2;[Pxx,f]=plomb(double(L1),Fs,Fmax,10,'power');figure(7);plot(Pxx);grid [Pxx_pks,f0]=findpeaks(Pxx,f,'MinPeakHeight',400)

depending upon the line chosen, when placing marker on spike tip, just above 300, there are 2 marker positions on top of the spike, there is need for higher frequency resolution. You choose a different line and get a better result, just one sample on tip of spike achieved without telling plomb() what Fs and Fmax to use:

[Pxx,f]=plomb(double(L1),nL1);figure(7);plot(Pxx);grid on

marker on f=123, started all over with different line, now length(L1)=379 tried with different windows, Welch for instance broadens peak of interest, and the frequency spike should be as sharp as possible insert071 insert072

again count cycles:

123/758*379/2 % =30.75

Nyquist wrote: 9/10 you try something, there is a better way:

14.- correlate signal with expected pulse, have a look at a single received pulse

try different pulses, gradually make the pulse narrower and taller

pulse1=[0 200 200 200 200 200 200 200 200 200 200 200 0] % poor correlation
pulse2=[0 50 100 100 100 100 100 100 100 100 100 50 0] % a bit better
pulse3=[0 0 50 100 100 100 100 100 100 100 50 0 0] % a bit better, just missing 1 pulse
pulse4=[0 0 0 50 100 100 100 100 100 50 0 0 0]
pulse5=[0 0 0 100 200 200 200 200 200 100 0 0 0]
pulse6=[0 0 0 0 100 200 200 200 100 0 0 0 0] % the count of peaks is now correct, do pulse1=[0 0 0 0 200 400 400 400 200 0 0 0 0]

pulse_test=[pulse1;pulse2;pulse3;pulse4;pulse5;pulse6]

for k=1:1:6
  r_L1_pulse1=conv(L1,pulse_test(k,:))
  hold on;figure(5);plot(r_L1_pulse1)
end
grid on

if you manually check with findpeaks, at every pulse there is improvement. Yet, depending on chosen sample line, there are missed pulses. If you start the sample line too early and let it run too late you may have to remove the first and last peak counts

    prob2=[]
    for k=1:20
   [pks_xcorr1puls,locs_xcorr1puls]=findpeaks(r_L1_pulse1,'minpeakdistance',k)
      prob2=[prob2 numel(pks_xcorr1puls)]
      end
      % figure(9);plot(prob2);grid on
     [pks_xcorr1puls,locs_xcorr1puls]=findpeaks(r_L1_pulse1,'minpeakdistance',6)
      % remove start and end counts
    numel(pks_xcorr1puls)‐2 % ‐2, as mentioned, for first and last of previous picture
    % 1st half of r_L1_pulse1 does not have relevant information
    r_L1_pulse1=r_L1_pulse1([375:750])
    n_r_L1_pulse1=[1:1:length(r_L1_pulse1)]

[Pxcorr,f]=plomb(r_L1_pulse1,n_r_L1_pulse1);figure(10);plot(Pxcorr);grid on
  length(Pxcorr)=752
122/752*376/2 % =30.5 pretty close, but still not as accurate as counting peaks

spectrum of the autocorrelation

15.- Improve by using command FILTER input signal before XCORR. Use the wonderful FDATOOL to build a direct FIR band pass filter equiripple that allows through the ribbons frequency only. Try different normalized BPF stop1 pass stop2 configurations and use 'a' 'b' to filter the sample line b: numerator a: denominator Mind the gap, once came across a book using b:den a:num

I found normalized band pass filter shoulders [.25 .3] and [.4 .45] reject 30dB enough Tried to attach ribbon_filter.fda for the readers to import it directly into FDATOOL, but the browser editor refused arguing it was none of the file formats accepted, and .fda files are not accepted as attachments (!?).

Could some one in Mathworks please make it possible for .fda files to be attached to answers/comments in the MATLAB CENTRAL forum? thanks in advance.

So back on track, manually doesn't take that long once you know where are the band pass shoulders. Filter mask

 filter pole/zero diagram

correction: the previous mask is a BPF that does not correspond to the following b and a. It's just that the following LPF is the first filter I got hands on in fdatool, it did the job, and did not update the coefficients, and now, moths latter after publishing this answer I realized. In any case if a LPF hits the mask, at least on nuts and bolts it used to be cheaper than BPF:

b=[0.028 0.053 0.071 0.053 0.028]
a=[1.000 ‐2.026 2.148 ‐1.159 0.279]
L1_BPfilt=filter(b,a,L1)
nL1_BPfilt=[1:length(L1_BPfilt)]
hold all;figure(1);plot(nL1,L1,nL1_BPfilt,L1_BPfilt);grid on
Comparing sample line without (blue) and with (orange) filtering:

16.- Improve by using command DETREND: removes VLF. Detrend the sample line before filtering

    L1_detrended=detrend(L1,'linear',[95 289]);hold all;figure(1);plot(L1_detrended)

detrend before filtering

detrend after filtering

L1_BPfilt_detr=detrend(L1_BPfilt,'linear',[95 289])

now with

plot(abs(fft(L1_BPfilt_detr)))

the highest spikes are around the area of interest, the frequency of the ribbons.

Yet depending on the sample line chosen, the top spike may be the right one, or you may erroneously choose a nearby taller spike by mistake.

length(abs(fft(L1_BPfilt_detr))) % =379

since we don't know the amount of ribbons in advance we cannot make the filter too narrow otherwise we might reject the right frequency.

17.- Improve by autocorrelating all received signal on itself, not on just a single pulse.

r_L1=xcorr(L1_BPfilt_detr,L1_BPfilt_detr)
r_L1=r_L1([374:end]) % self correlations are symmetric, we only one one side, taking right figure(2);plot(r_L1);grid on

haven't done it here, I amplify a bit 3*(G+B), but the recommended place to amplify is AFTER filtering, not before. After checking different FINDPEAKS 'minpeakheight', all and each ribbon has 1 and only 1 peak above 0 when 'minpeakheight'=0

[pks,locs]=findpeaks(r_L1,'minpeakheight',0)
numel(pks) % = 29

note that this graph still includes start stop as peaks, 312

18.- So far it is only 1 line, let's remove luck out of the loop, with a for loop repeat for all lines in the sample strip. Recap:

ribbon_count=zeros(1,40)
for k=400:440 
im_line_2=A(k,:,2);im_line_3=A(k,:,3) 
im_line_2=255‐im_line_2;im_line_3=255‐im_line_3 % invert 
L1=3*(double(im_line_2)+double(im_line_3)) % 3*(L1=L1‐mean(L1) 
L1_BPfilt=filter(b,a,L1) 
L1_BPfilt_detr=detrend(L1_BPfilt,'linear',[95 289]) 
r_L1=xcorr(L1_BPfilt_detr,L1_BPfilt_detr);r_L1=r_L1([374:end]) % XCORR: 
[pks,locs]=findpeaks(r_L1,'minpeakheight',0) 
ribbon_count(k‐400+1)=numel(pks) 
end

have a look

stem(ribbon_count)
hist(ribbon_count)

the answer is

mode(ribbon_count) % overwhelming majority to: 29

If you find this answer of any help solving the question above, please click on the thumbs-up vote link, thanks in advance

John

[e-mail]

UPC-ETSETB 99

comment 1: there are other tools and ways, for instance medfilt1, sgolay, sgolayfilt: all seem to be able to remove HF. I was about to try sgolay, and medfilt1, that both seem to remove HF, but fdatool went well enough calculating 'a' and 'b'.

comment 2: you see that some pulses have 'horns' and many have slope on top, and between pulses there are 'little pulses'? that can be caused by multipath, aka echo.

Within certain limitations equalizers remove echo.

But equalizers need training sequences that precisely should have the right frequency, which is what the question is after

There is a built in class to equalize signals.

I wanted to use the BER calculator with the QAM carrier sync recovery, by sweeping a frequency range with:

• guess a frequency
• build training sequence
• clean all lines in strip with equalizer
• filter, amplify, ..
• check BER

It would be reasonable to expect a BER notch on the right amount of ribbons but it would have taken far longer that the above lines.

[John BG]

so, the popular saying, if a LPF does the job no one will care to replace it with a BPF, it's true even in MATLAB ;)

[Walter Roberson]

"Could some one in Mathworks please make it possible for .fda files to be attached to answers/comments in the MATLAB CENTRAL forum?"

You can zip the files and attach the .zip

[John BG]

Thanks Walter

Answer 3

[Selva Karna]

Thank you so much sir.....

[John BG]

Hi Selva

I had already put this one in the cold cases folder.

The green board with the parallel strips is a prototype of antenna, isn't it?