Construire une image RVB (à partir image multispectrale)

**jean-pat** · 23/05/2011, 11h16

Bonjour à tous,
J'essaie de construire une image RVB à partir de trois images monochromes (des tableaux numpy.array):

R, G, B

. J'ai suivi les indications de l'auteur de astrobetter.Le code est censé produire un tableau "rgb".
Dans mon code cela donne:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#Reconstitute original image size
#should be a red image of size =shape
R=c3_r.reshape(shape)
G=c3_g.reshape(shape)
B=c3_b.reshape(shape)
#Try to build a rgb image
rgb = np.zeros((shape[0],shape[1],3),dtype=float)
print rgb.shape
mxr=np.max(R)
mxg=np.max(G)
mxb=np.max(B)
#normalize to 0-255 uint8
Rnorm=np.uint8((255*(R/mxr)))
Gnorm=np.uint8((255*(R/mxg)))
Bnorm=np.uint8((255*(R/mxb)))
#copy each RGB component in an RGB array
rgb[:,:,0]=Rnorm
rgb[:,:,1]=Gnorm
rgb[:,:,2]=Bnorm

Quand j'essaie d'afficher cette image avec:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
pylab.gray()
pylab.subplot(224,frameon=False, xticks=[], yticks=[])
pylab.imshow(rgb)
pylab.show()

Je ne vois pas une image rvb mais une image en niveaux de gris avec un fond blanc (au lieu de noir) et je n'arrive pas à la sauver en png.

Pour raconter toute l'histoire, j'ai au départ cinq images correspondant une scène vue à différentes longueur d'ondes (du bleu au proche infra rouge). Je souhaite combiner ces cinq images en trois pour les visualiser simultannement. On peut faire cela en quelques clic avec imagej et le plugins image5D.
J'essaie de refaire la même chose en python+numpy+pylab. J'aurais également des questions sur la construction de la matrice m.Voila le code complet:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
import numpy as np
import readmagick
import os
import pylab
user=os.path.expanduser("~")
workdir=os.path.join(user,"Applications","ImagesTest","MFISH")
blue="Aqua.tif"
green="Green.tif"
gold="Gold.tif"
red="Red.tif"
frd="FarRed.tif"
complete_path=os.path.join(workdir,blue)
i1=readmagick.readimg(complete_path)
#
complete_path=os.path.join(workdir,green)
i2=readmagick.readimg(complete_path)
#
complete_path=os.path.join(workdir,gold)
i3=readmagick.readimg(complete_path)
#
complete_path=os.path.join(workdir,red)
i4=readmagick.readimg(complete_path)
#
complete_path=os.path.join(workdir,frd)
i5=readmagick.readimg(complete_path)
#
shape=i5.shape
 
b=i1.flatten()
g=i2.flatten()
y=i3.flatten()
r=i4.flatten()
f=i5.flatten()
#make a 2D array
#each line is a pixel, column=channel
channel5=np.vstack((f,r,y,g,b))
c5=channel5.T
#print c5.shape
#try some coef
m=np.array([[0.80, 0.10, 0.01],
            [0.10, 0.80, 0.04],
            [0.10, 0.05, 0.10],
            [0.00, 0.03, 0.15],
            [0.00, 0.02, 0.70]])
print m.shape
 
#produce  a flat rgb image
#combine linearly the five channels into three RGB channels
c3=np.dot(c5,m)
 
#isolate R,G,B component
c3_r=c3[:,0]
c3_g=c3[:,1]
c3_b=c3[:,2]
#Reconstitute original image size
#should be a red image of size =shape
R=c3_r.reshape(shape)
G=c3_g.reshape(shape)
B=c3_b.reshape(shape)
#Try to build a rgb image
rgb = np.zeros((shape[0],shape[1],3),dtype=float)
print rgb.shape
mxr=np.max(R)
mxg=np.max(G)
mxb=np.max(B)
#normalize to 0-255 uint8
Rnorm=np.uint8((255*(R/mxr)))
Gnorm=np.uint8((255*(R/mxg)))
Bnorm=np.uint8((255*(R/mxb)))
#copy each RGB component in an RGB array
rgb[:,:,0]=Rnorm
rgb[:,:,1]=Gnorm
rgb[:,:,2]=Bnorm
#complete_path=os.path.join(workdir,"R.png")
#readmagick.writeimg(R,complete_path)
pylab.gray()
#
pylab.subplot(221,frameon=False, xticks=[], yticks=[])
pylab.imshow(Rnorm)
#complete_path=os.path.join(workdir,"R.png")
#readmagick.writeimg(Rnorm,complete_path)
#
pylab.subplot(222,frameon=False, xticks=[], yticks=[])
pylab.imshow(Gnorm)
#complete_path=os.path.join(workdir,"G.png")
#readmagick.writeimg(Gnorm,complete_path)
#
pylab.subplot(223,frameon=False, xticks=[], yticks=[])
pylab.imshow(Bnorm)
#complete_path=os.path.join(workdir,"B.png")
#readmagick.writeimg(Bnorm,complete_path)
#pylab.clf()
pylab.subplot(224, aspect='equal',frameon=False, xticks=[], yticks=[])
pylab.imshow(rgb)
pylab.show()

**jean-pat** · 23/05/2011, 19h34

J'ai simplifié mon code vite fait mal fait sans algèbre lineaire.

Le résultat n'est surrement pas optimal

**Alexis.M** · 22/06/2011, 10h29

Bonjour,

je pars du principe que tu as récupéré tes images dans i1,i2,...,i5:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
multi_spectr = np.dstack([i1,i2,i3,i4,i5])
conv_matrix = np.array([
...
])
 
shape = multi_spectr.shape
multi_spectr = multi_spectr.reshape(shape[0]*shape[1],shape[2])
rgb = np.dot(multi_spectr,conv_matrix.T).reshape(shape[:2]+(3,))

**jean-pat** · 22/06/2011, 13h35

merci, je vais essayer cela.

edit:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
# -*- coding: utf-8 -*-
"""
Created on Mon May 23 14:07:40 2011
 
@author: jean-pat
"""
 
import numpy as np
import os
import pylab
#from scipy import linalg as la
#from scipy import ndimage as nd
import scipy as sp
user=os.path.expanduser("~")
workdir=os.path.join(user,"Applications","ImagesTest","MFISH")
blue="Aqua.tif"
green="Green.tif"
gold="Gold.tif"
red="Red.tif"
frd="FarRed.tif"
complete_path=os.path.join(workdir,blue)
i1=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,green)
i2=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,gold)
i3=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,red)
i4=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,frd)
i5=sp.misc.imread(complete_path)
######################
multi_spectr = np.dstack([i1,i2,i3,i4,i5])
conv_matrix = np.array([
[0.8,0.15,0.05,0,0],
[0,0.05,0.05,0.8,0.1],
[0,0,0.05,0.15,0.8]
])
 
shape = multi_spectr.shape
multi_spectr = multi_spectr.reshape(shape[0]*shape[1],shape[2])
rgb = np.dot(multi_spectr,conv_matrix.T).reshape(shape[:2]+(3,))
rgb8=np.uint8(rgb)
########################
#R=0.8*i5+0.15*i4+0.05*i3+0.00*i2+0.00*i1
#V=0.00*i5+0.05*i4+0.05*i3+0.80*i2+0.10*i1
#B=0.00*i5+0.00*i4+0.05*i3+0.15*i2+0.80*i1
#shape=i5.shape
#rgb = np.zeros((shape[0],shape[1],3),dtype=np.uint8)
#rgb[:,:,0]=np.copy(R)
#rgb[:,:,1]=np.copy(V)
#rgb[:,:,2]=np.copy(B)
pylab.subplot(111,frameon=False, xticks=[], yticks=[])
pylab.imshow(rgb8)
pylab.show()

Il faut convertir en 8 bits pour voir l'image avec pylab.

merci

**jean-pat** · 22/06/2011, 14h10

Au passage, existe t-il de "meilleures" matrices conv_matrix que d'autres?
A mon sens une bonne matrice est une matrice qui permet au mieux de séparer les couleurs.
On peut imaginer qu'il en existe une bonne qui donne une image RVB, cette image peut ensuite être convertie en HSI. Par rotation de la composante hue on peut d'une image à une autre. Il devrait donc y avoir "beaucoup" de matrices "équivalentes" parce qu'il est possible de passer d'une image RGB à l'autre en modifiant la composante Hue.

Initialement, les cinqs images de départ codent pour 24 couleurs différentes. Peut pratiquement faire une ACP (SVD...) pour construire cette image RVB?

merci

Jean-Patrick

**Alexis.M** · 22/06/2011, 14h34

Effectivement une ACP aurait l'avantage d'être optimale pour tes données, tu peux imaginer prendre une ou plusieurs photos pour calculer la matrice de conversion une fois pour toute.

**jean-pat** · 22/06/2011, 18h00

malheureusement, il ne s'agit pas d'images que j'ai prise, je ne dispose pas du temps d'exposition de chaque image. Il s'agit de données publiées il y a déjà du temps:

Copyright (c) 2000 Advanced Digital Imaging Research
All rights reserved.

Permission is hereby granted, without written agreement and without license or royalty fees, to use, copy, modify, and distribute these images and their documentation for any purpose, provided that the above copyright notice and the following two paragraphs appear with all copies of these images.

IN NO EVENT SHALL PSI RESEARCH BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF THESE IMAGES AND THEIR DOCUMENTATION, EVEN IF PSI RESEARCH HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

PSI RESEARCH SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE IMAGES PROVIDED HEREUNDER ARE ON AN "AS IS" BASIS, AND PSI RESEARCH HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.

--------------------------------------------------------------------------------

PSI M-FISH IMAGE SET #1

The accompanying image set contains 200 M-FISH Images. The images are divided into directories based on karyotype and slide.

File Nomenclature

The first character in the directory name represents the probe set (A, ASI; P, PSI/Cytocell; V, Vysis). The second two characters are a number designating which slide the images in that directory came from. The next two characters are underscores, and the final two characters are a code for the karyotype of that image (The ISCN designations for each code are listed below.) Therefore a directory containing only the images with normal male karyotypes from slide 99 (using PSI probes) would be P99_XY.

The names of the images follow a similar naming scheme. The first character represents the probe set. The next two characters represent the slide number that the image came from. The next two characters are the number of that image on the slide. The final two characters represent the karyotype code. Therefore, if image number 12 from slide 98 (using ASI probes) were a normal female metaphase, its file name would be A9812XX.

M-FISH Format

If the image set is stored in PSI's M-FISH format, all probes and image planes are stored in the same file. If a file ends in an E, it means that that image has been enhanced manually. Any file that ends in E, unless it is an extreme case (labeled karyotype code EX), also contains my diagnosis of the correct karyotype for that metaphase.

If the images set is stored in PSI's M-FISH format, one directory from each slide will also contain the M-FISH classifier used to come up with the karyotype.

Non-MFISH Format (PNG, TIFF, JPEG, PICT)

If the images are not stored in PSI's M-FISH format, each plane of the image is stored in a separate file. The file name will have 8 characters, and the 8th and last character represents the probe that the image represents.

3 532
5 Cy 5.5
6 568
A S. Aqua
C Cy 5
D DAPI
E DEAC
F Far Red
G S. Green
I FITC
O S. Orange
R S. Red
T Texas Red
Y S. Gold

In addition to a file for each plane in the M-FISH image, there is another file that contains my diagnosis of the correct karyotype for that metaphase. This file has a 'K' for its 8th character. This image is labeled such that all the pixels belonging to chromosome 1 are 1, all the pixels belonging to chromosome 2 are 2, etc. All background pixels are 0, and all pixels that are part of a region of overlap are -1. (This is a little meaningless in the case of a translocation. In this case, the whole chromosome is labeled the same-- generally that class which makes up the most of the chromosome.)

In non-MFISH versions of this image set, each directory will also have a text file (of the same name as the directory) which contains the ISCN designation of the karyotype contained in that directory.

DISCLAIMER:

The karyotypes in this image set are my interpretation of the metaphases. If you feel that I have made an error in my diagnosis of any image, please let me know at wade@psires.com

August 25, 2000
Advanced Digital Imaging Research
Wade Schwartzkopf
wade@psires.com

**Alexis.M** · 22/06/2011, 18h54

De toute façon il n'est pas très difficile de tester en utilisant les données d'une ou plusieurs images. ça n'a pas moins de sens qu'utiliser une matrice de transformation arbitraire...

**jean-pat** · 22/06/2011, 22h46

Envoyé par Alexis.M

Bonjour,

je pars du principe que tu as récupéré tes images dans i1,i2,...,i5:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
multi_spectr = np.dstack([i1,i2,i3,i4,i5])
conv_matrix = np.array([
...
])
 
shape = multi_spectr.shape
multi_spectr = multi_spectr.reshape(shape[0]*shape[1],shape[2])
rgb = np.dot(multi_spectr,conv_matrix.T).reshape(shape[:2]+(3,))