Đồ án Nghiên cứu mô phỏng nguyên lý hoạt động của thiết bị chụp cắt lớp cộng hưởng từ hạt nhân

kết quả, cho phép lập trình với giao diện làm việc rất dễ dàng cho người sử dụng. Dữ liệu cùng với thư viện được lập trình sẵn cho phép người sử dụng có thể có được những ứng dụng sau đây: • Sử dụng các hàm có sẵn trong thư viện, các phép tính toán học thông thường • Cho phép lập trình tạo ra những ứng dụng mới. • Cho phép mô phỏng các mô hình thực tế. • Phân tích, khảo sát và hiển thị dữ liệu. • Có phần mềm đồ hoạ cực mạnh. • Cho phép phát triển, giao tiếp với một số ngôn ngữ lập trình khác như C++, Fortran. b. Tổng quan về cấu trúc dữ liệu của MATLAB, các ứng dụng Matlab là một hệ thống tương giao, các phần tử dữ liệu là một mảng (mảng này không đòi hỏi về kích thước). Chúng cho phép giải quyết các vấn đề liên 53 quan đến lập trình bằng máy tính, đặc biệt sử dụng các phép tính về ma trận hay véctor và có thể sử dụng ngôn ngữ C học Fortran lập trình rồi thực hiện ứng dụng lập trình đó bằng các câu lệnh goị từ MATLAB. MATLAB được viết tắt từ chữ MATRIX LABORATORY tức là thư viện về ma trận, từ đó phần mềm MATLAB được viết nhằm cung cấp cho việc truy cập vào phần mềm ma trận một cách dễ dàng, phần mềm ma trận này được phát triển bởi các công trình Linpack và Eispack. Ngày nay MATLAB được phát triển bởi Lapack và Artpack tạo nên một nghệ thuật phần mềm cho ma trận. *Dữ liệu Dữ liệu của MATLAB thể hiện dưới dạng ma trận (hoặc mảng –tổng quát), và có các kiểu dữ liệu được liệt kê sau đây: • Kiểu đơn single, kiểu này có lợi về bộ nhớ dữ liệu vì nó đòi hỏi ít byte nhớ hơn, kiểu dữ liệu này không được sử dụng trong các phép tính toán học, độ chính xác kém hơn. • Kiểu double kiểu này là kiểu thông dụng nhất của các biến trong MATLAB. • Kiểu Sparse. • Kiểu int8, uint8, int16... • Kiểu char ví dụ ‘Hello’ • Kiểu cell. • Kiểu Structure. Trong MATLAB kiểu dữ liệu double là kiểu mặc định sử dụng trong các phép tính số học. Các bạn có thể tham khảo các kiểu dữ liệu khác trong đĩa CD Help MATLAB 7.0 * ứng dụng MATLAB tạo điều kiện thuận lợi cho: 54 • Các khoá học về toán học • Các kỹ sư, các nhà nghiên cứu khoa học • Dùng MATLAB để tính toán, nghiên cứu tạo ra các sản phẩm tốt nhất trong sản xuất. • Giải quyết các thuật toán mô phỏng trong kỹ thuật một cách hiệu quả và chính xác. *Toolbox Ngoài ra MatLab còn cung cấp một hệ thống thư viện các hộp công cụ (Toolbox) cho phép bạn ứng dụng các kỹ thuật, lập trình để phân tích, thiết kế, mô phỏng các mô hình hay các bài toán kỹ thuật từ đơn giản đến phức tạp. Một số hộp công cụ có thể kể tên ở đây như: Simulink, xử lý tín hiệu, xử lý ảnh, logic mờ, mạng nơron,vv. Trở lại với mục tiêu đề tài của ta là: “Nghiên cứu mô phỏng nguyên lý hoạt động của thiết bị chụp cắt lớp CHTHN”. Từ việc nghiên cứu nguyên lý hoạt động của thiết bị này trên tất cả các khâu xử lý cả bằng phần cứng lẫn phần mềm trong chương 1 và chương 2, ta nhận thấy các tín hiệu xử lý trên các khâu phần cứng của thiết bị đều có thể biểu diễn dưới dạng vector, ma trận. Ngoài ra các khâu xử lý bằng phần mềm của thiết bị có thể được thực hiện bằng các thuật toán thông dụng như phép biến đổi Fourier, tích chập, các hàm lượng giác với các dữ liệu ở dạng vector, ma trận. Với các đặc điểm đã xét ở trên về ngôn ngữ lập trình MatLab và kết hợp với việc khảo sát thêm một số công cụ phần mềm thông dụng khác như: Visual C++, Visual Basic, Labview, Pascalcó thể kết luận việc lựa chọn ngôn ngữ lập trình Matlab để thực hiện việc mô phỏng là thích hợp hơn cả. Tuy nhiên với tất cả những phần mềm này đều có điểm hạn chế là với các tín hiệu có tần số rất cao (ở đây tín hiệu CHTHN có dải tần từ vài chục MHz đến vài trăm MHz) thì việc lấy mẫu tín hiệu và xử lý rất khó khăn. Hạn chế này sẽ 55 được xem xét kỹ hơn trong phần mô phỏng hoạt động của bộ tách sóng pha cầu phương sau này. 56 CHƯƠNG 3 MÔ HÌNH TOÁN HỌC VÀ SƠ ĐỒ THUẬT TOÁN MÔ PHỎNG NGUYÊN LÝ HOẠT ĐỘNG CỦA THIẾT BỊ CHỤP CẮT LỚP CỘNG HƯỞNG TỪ HẠT NHÂN Như đã phân tích trong chương 2, nguyên lý hoạt động của thiết bị CCL CHTHN được thể hiện ở từng công đoạn xử lý thông tin bằng cả phần cứng lẫn phần mềm. Vì vậy việc mô phỏng cần được tiến hành lần lượt với từng phần (từng công đoạn xử lý tín hiệu) rồi sau đó sẽ kết hợp các phần lại để hình thành chương trình mô phỏng cho toàn bộ thiết bị. 3.1. Mô phỏng từng công đoạn xử lý của thiết bị Dựa trên những nghiên cứu tổng quan về nguyên lý hoạt động của cả phần cứng lẫn phần mềm của thiết bị chụp cắt lớp CHTHN trong các chương 1 và 2, có thể phân tích toàn bộ hệ thống thành các khâu xử lý tín hiệu CHTHN nối tiếp và tương ứng với chúng là các công đoạn mô phỏng như sau: ●Mô phỏng quá trình kích thích và tạo giả tín hiệu CHTHN. ●Mô phỏng quá trình xử lý tín hiệu CHTHN bằng bộ tách sóng pha cầu phương. ●Mô phỏng quá trình xử lý dữ liệu và hiển thị ảnh bằng phần mềm. Trong công đoạn đầu tiên ta thực hiện mô phỏng quá trình kích thích và tạo giả tín hiệu CHTHN phát ra từ các kiểu mô khác nhau ở một cường độ từ trường cụ thể. Tiếp theo các tín hiệu CHTHN giả sẽ được đưa đến xử lý bởi bộ tách sóng pha cầu phương. Tại đây tín hiệu vô tuyến tần số cao (hàng trăm MHz) sẽ được chuyển sang dạng tín hiệu âm tần có tần số thấp. Tín hiệu thu được sau khâu này được tách thành hai thành phần, thường được gọi là phần “thực” và phần “ảo” tương ứng với hai thành phần từ trường vuông góc nhau và . Bằng cách sử dụng bộ tách sóng cầu phương có hai kênh “thực” XM YM 57 và “ảo” như vậy cho phép dễ dàng xác định được biên độ và pha của mỗi tín hiệu. Ở công đoạn mô phỏng cuối cùng sẽ tiến hành xử lý dữ liệu thu được sau bộ tách sóng pha cầu phương bằng thuật toán biến đổi Fourier 2 chiều và thực hiện các bước hiệu chỉnh cần thiết để nhận được ảnh CCL CHTHN mong muốn. 3.1.1. Mô phỏng quá trình kích thích và tạo giả tín hiệu CHTHN Chúng ta biết rằng sự kích thích và chương động của vectơr từ trường tổng hợp tuân theo phương trình Bloch (1.9). Lời giải của nó đã được biết đến trong phần lý thuyết về hiện tượng CHTHN bao gồm sự tổng hợp của các quá trình dãn hồi T1, T2 và quá trình chương động tự do. Như vậy thực hiện mô phỏng phương trình Bloch là cách chính xác nhất để nghiên cứu tác động kích thích của một dãy xung vô tuyến lên vector từ trường. Thực chất của quá trình chương động tự do (hay chính là sự hình thành tín hiệu CHTHN) là việc quay vector từ trường quanh trục Z kết hợp với sự suy giảm của nó theo các hằng số dãn hồi T1,T2. → M Z Y X Hình 3.1: Quá trình chương động tự do Trước khi đưa ra mô hình toán học và sơ đồ thuật toán mô phỏng ta cần có một số quy ước sau: 1. Biểu diễn từ trường dưới dạng ma trận 3x1. 2. Biểu diễn các thời gian dãn hồi T1, T2 dưới dạng ma trận. 58 3. Dùng các hàm để thực hiện việc quay vector từ trường. 4. Xác định kết quả của sự dãn hồi tự do (sự chương động) tổng hợp. Theo kết quả giải phương trình Bloch ở chương 1 các nghiệm được xác định theo công thức (1.13) : ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ −−= −= −= ))T texp((1M(t)M t)sin(ω)MT texp((t)M t)cos(ω)MT texp((t)M 1 0z o0 2 y o0 2 x (3.1) Ở đây )0()0( 220 yx MMM += Khi này ta có thể biểu diễn : ●Quá trình suy giảm của vector từ trường tổng hợp theo cả hai hằng số thời gian T1 và T2 là: Mdh=A*M(0)+B (3.2) Ở đây Mdh là vector biểu diễn sự suy giảm của vector từ trường tổng hợp. M(0) là vector biểu diễn giá trị ban đầu của từ trường tổng hợp kể từ thời điểm kích thích. Giả thiết M(0)=[M0 0 0]’ A và B là các ma trận, vector biểu diễn cho việc kích thích và quá trình chương động tự do của từ trường. A=[exp(-t/T1), 0, 0;0, exp(-t/T1), 0; 0, 0, exp(-t/T2)] (3.3) B=[0, 0, M0-M0exp(-t/T2)]' (3.4) Dấu “ ' ” biểu thị dạng ma trận chuyển vị. ●Quá trình chương động (quay) quanh trục Z của vector từ trường tổng hợp tuân theo: Mrot = RotZ.M (3.5) Ở đây Mrot là vector biểu diễn cho quá trình quay vector từ hoá quanh trục Z. 59 M là vector từ hoá tổng hợp. RotZ là ma trận thực hiện việc quay M r quanh trục Z. RotZ(ϕ ) = [ ϕcos - ϕsin 0 ; ϕsin ϕcos 0 ; 0 1] (3.6) 0 Khi thực hiện kích thích vô tuyến ta có thể quay vector từ trường M r theo các phương trục X, Y bằng các ma trận sau: RotX(ϕ ) = [1 0 ; 0 0 ϕcos - ϕsin ; 0 ϕsin ϕcos ] (3.7) RotY (ϕ )=[ ϕcos -0 ϕsin ; 1 ; -0 0 ϕsin 0 ϕcos ] (3.8) Tổng hợp của hai quá trình này ta có quá trình chương động tự do của vector từ trường: Mcdtd = Acd*M(0)+Bcd (3.9) Ở đây: Mcdtd là vector biểu diễn cho quá trình chương động tự do của vector từ hoá. M(0) là giá trị ban đầu của vector từ hoá. Với M(0)=[M0 0 0]’. Acd, Bcd là các vector, ma trận thực hiện quá trình chương động. Acd = [exp(-t/T2) 0 0; 0 exp(-t/T2) 0; 0 0 exp(-t/T1)]*RotZ(ϕ ) (3.10) Bcb = [0 0 M0-M0exp(-t/T1)]’ (3.11) Góc ϕ xác định bởi tần số xung vô tuyến kích thích theo công thức (1.27): τπϕ ..2 f= ( là tần số xung vô tuyến, f τ là khoảng thời gian kích thích). Với từng phương pháp tạo ảnh cụ thể thì ta có các dãy xung kích thích khác nhau. Ở đây chúng ta sẽ nghiên cứu một vài phương pháp tạo ảnh thường được dùng như: Phương pháp tiếng vọng spin, phương pháp phục hồi đảo ngược, và phương pháp tiếng vọng gradient. ● Phương pháp tiếng vọng spin 60 Hình 3.2: Biểu đồ dãy xung kích thích tiếng vọng spin Trong phương pháp này với mỗi chu kỳ thu nhận ta thực hiện hai kích thích. Đầu tiên thực hiện quay vector từ trường xuống mặt phẳng ngang XY một góc 900 so với trục Z, khi này bắt đầu xẩy ra quá trình chương động tự do của các thành phần từ trường. Sau một khoảng thời gian TE/2 thực hiện quay tiếp đi một góc 1800 và tại thời điểm TE sẽ nhận được tín hiệu tiếng vọng. ●Phương pháp phục hồi đảo ngược. Hình 3.3: Biểu đồ dãy xung kích thích phục hồi đảo ngược Phương pháp tạo ảnh này khác biệt so với phương pháp tiếng vọng spin ở chỗ tại đầu mỗi chu kỳ thu nhận TR ta sẽ dùng một xung đảo ngược quay vector từ trường đi một góc 1800, rồi sau một khoảng thời gian TI sẽ thực hiện các kích thích như đối với phương pháp tiếng vọng spin để thu tín hiệu tiếng vọng. Khi này quá trình dãn hồi của vector từ hoá bắt đầu từ giá trị âm chứ không phải từ giá trị 0 như trong phương pháp tiếng vọng spin. ●Phương pháp tiếng vọng Gradient 61 Hình 3.4: Biểu đồ dãy xung kích thích tiếng vọng gradient Phương pháp tiếng vọng gradient cũng giống như phương pháp tiếng vọng spin nhưng lúc này xung kích thích ban đầu tạo góc quay nhỏ hơn 900 và chu kỳ thu nhận TR thường rất ngắn. Vì vậy phương pháp này còn được gọi là phương pháp kích thích với góc lệch nhỏ. Trên cơ sở mô hình toán học của công đoạn kích thích và tạo tín hiệu CHTHN có thể xây dựng sơ đồ thuật toán chung mô phỏng được chỉ ra trong hình 3.5: Quá trình mô phỏng diễn ra như sau: Trên sơ đồ thuật toán mỗi tín hiệu tương ứng với một kiểu mô cụ thể sẽ có các tham số cần thiết như: Thời gian dãn hồi T1, T2, mật độ proton. Ở cường độ từ trường đủ lớn đảm bảo phát sinh hiện tượng CHTHN, tiến hành kích thích bằng cách quay vector đi các góc cần thiết (Với từng phương pháp tạo ảnh sẽ có các kiểu kích thích tương ứng) sau khi kết thúc sẽ xẩy ra quá trình dãn hồi tự do. Thực hiện lưu giữ các giá trị thay đổi theo thời gian của cả 3 thành phần từ trường , , . Dữ liệu này biểu diễn dạng tín hiệu CHTHN, và sau đó được xử lý, hiển thị bằng một số hàm đồ hoạ trong Matlab. XM YM ZM 62 Hình 3.5: Sơ đồ thuật toán mô phỏng quá trình kích thích và tạo tín hiệu CHTHN. 3.1.2. Mô phỏng quá trình xử lý tín hiệu bằng bộ tách sóng pha cầu phương Tách sóng pha cầu phương là một phương pháp giải điều chế tín hiệu thông dụng không chỉ được dùng trong thiết bị CCL CHTHN mà còn được dùng trong nhiều lĩnh vực khoa học kỹ thuật khác. 63 Trong CCL CHTHN, kỹ thuật này đã được nói đến trong phần 1.3 của chương 1. Về cơ bản hoạt động của nó có thể được mô tả theo mô hình toán học sau: Giả sử tín hiệu đi vào bộ tách sóng cầu phương là tín hiệu dạng sin có tần số ω , tín hiệu chuẩn cũng là tín hiệu dạng sin nhưng có tần số 0ω . Khi đó tín hiệu đầu vào bộ tách sóng cầu phương có dạng : )2cos( tπω Và tín hiệu chuẩn đưa vào 2 bộ trộn có dạng: )2cos( 0tπω và )2sin( 0tπω Theo công thức biến đổi lượng giác, đầu ra hai bộ trộn là : )2cos( tπω * )2cos( 0tπω = 1/2 )22cos( 0tt πωπω − + 1/2 )22cos( 0tt πωπω + )2cos( tπω * )2sin( 0tπω = 1/2 )22sin( 0tt πωπω − + 1/2 )22sin( 0tt πωπω + Khi đi qua bộ lọc thông thấp thành phần tần số tổng bị loại bỏ, chỉ còn lại các thành phần tần số hiệu: Đầu ra 1 = 1/2 )22cos( 0tt πωπω − Đầu ra 2 = 1/2 )22sin( 0tt πωπω − Trong thiết bị CCL CHTHN tín hiệu đưa đến bộ tách sóng pha cầu phương là tín hiệu cộng hưởng từ được cảm ứng trong các cuộn vô tuyến RF. Theo công thức (1.29) dạng tín hiệu CHTHN đưa vào có dạng: Ở đây tín hiệu S(t) là tổng của tất cả các voxel phát ra từ một lát cắt đã được chọn sẵn. )(rω , )(rϕ là các giá trị tần số và pha đã được mã hoá theo vị trí. Nếu tín hiệu chuẩn có dạng )2cos( tRπω thì tín hiệu trên hai đầu ra bộ tách sóng sẽ là: 64 Ở đây độ lệch tần số : Rr ωωω −=∆ )( (3.15) Việc thu nhận tín hiệu CHTHN bằng bộ tách sóng pha cầu phương có vai trò đặc biệt quan trọng trong việc tái tạo các ảnh chụp cắt lớp ở các quá trình xử lý sau này. Đây là quá trình xử lý hoàn toàn bằng phần cứng và thực hiện với các tín hiệu CHTHN dạng tương tự. Điều này có thể được giải thích bởi đặc điểm của tín hiệu CHTHN là dạng tín hiệu vô tuyến có tần số rất cao (đến hàng trăm MHz). Nếu xử lý bằng phần mềm thì phải lấy mẫu để số hoá tín hiệu và thực hiện các quá trình xử lý số. Điều này yêu cầu phải có một hệ thống máy tính rất mạnh mới có thể đáp ứng được. Tuy nhiên đây chỉ là trên lý thuyết vì thực tế khi lấy mẫu tín hiệu có tần số cao như tín hiệu CHTHN là rất khó khăn vì các bộ chuyển đổi tương tự số không đáp ứng được tần số xung nhịp cao như vậy. Theo lý thuyết lấy mẫu của Nyquist, tần số lấy mẫu ít nhất phải lớn hơn hai lần tần số của tín hiệu. Nếu không thoả mãn điều kiện này thì sẽ xẩy ra hiện tượng nhiễu loạn hay còn gọi “Aliasing” như đã được nói đến trong chương 1. Vì việc mô phỏng hiện tượng CHTHN cũng như các công đoạn xử lý tín hiệu của thiết bị trên cả phần cứng lẫn phần mềm đều được thực hiện trên máy tính nên mọi tín hiệu đều phải được lấy mẫu và xử lý dưới dạng số. Để mô phỏng được chính xác các quá trình này thì tần số lấy mẫu cần thiết sẽ rất lớn, điều này đòi hỏi phải có các công cụ mô phỏng cho phép thực hiện được trong dải tần như vậy. Với điều kiện chương trình mô phỏng thực hiện trên các máy tính thông thường và công cụ mô phỏng là phần mềm toán học Matlab thì không thể đáp ứng được đúng theo yêu cầu thực tế. Tuy nhiên do 65 mục đích chủ yếu là mô phỏng nguyên lý hoạt động của thiết bị nên có thể tiến hành trong dải tần thấp hơn mà vẫn phản ánh đúng được các quá trình xử lý trong các khâu của thiết bị. Cụ thể ở đây các tín hiệu giả CHTHN được chọn có dải tần khoảng vài chục KHz và tần số lấy mẫu khoảng vài trăm KHz. Trong quá trình mô phỏng ta có thể nâng dải tần này lên cao hơn nhưng khi đó do hạn chế bởi các lý do đã nêu ở trên nên quá trình thực hiện sẽ rất lâu và có thể làm treo hệ thống. Trên cơ sở mô hình toán học và các phân tích đưa ra ở trên có thể xây dựng sơ đồ thuật toán mô phỏng cho công đoạn này được chỉ trong hình 3.7. Quá trình mô phỏng diễn ra như sau: Sơ đồ thuật toán mô phỏng cho phép đưa vào một hoặc nhiều tín hiệu CHTHN. Mỗi tín hiệu sẽ được lưu giữ dưới dạng vector; sau khi thực hiện quét hết các tín hiệu sẽ nhận được các vector tín hiệu. Tính tổng của các vector này cho ta vector biểu diễn tín hiệu tổng hợp của tất cả các tín hiệu đưa vào. Đây là tín hiệu vô tuyến có tần số cao, là tổng hợp của nhiều tín hiệu có các thành phần tần số và pha khác nhau. Tín hiệu này được đưa đến trộn với hai tín hiệu chuẩn lệch pha nhau 900 trên hai kênh riêng biệt. Tín hiệu đầu ra trên hai kênh này được đưa qua bộ lọc thấp tần và lưu giữ dưới dạng vector phần thực (Re) và phần ảo (Im). Lúc này có thể xác định được biên độ và pha của tín hiệu dựa vào các giá trị phần thực và phần ảo. Hai thành phần tín hiệu này được tổng hợp lại theo dạng phức và được gọi là tín hiệu suy giảm cảm ứng tự do (hay còn được gọi là tín hiệu FID). Đây là tín hiệu thấp tần có tần số nằm trong dải âm tần, tín hiệu này sẽ được lưu trữ lại dưới dạng vector và sau đó bằng các quá trình xử lý cần thiết để có thể quan sát được dưới nhiều dạng khác nhau. 66 Hình 3.7: Sơ đồ thuật toán mô phỏng quá trình xử lý tín hiệu bằng bộ tách sóng pha cầu phương. 67 3.1.3. Mô phỏng quá trình xử lý và hiển thị ảnh Sau khi thu nhận và lưu trữ được ma trận dữ liệu thô trên các khâu xử lý bằng phần cứng của thiết bị, các quá trình xử lý tiếp theo được thực hiện bằng phần mềm nhờ hệ thống máy tính chuyên dụng. Trong công đoạn này ta thực hiện mô phỏng các quá trình xử lý dữ liệu như: Biến đổi Fourier 2 chiều, các phép nội suy 2 chiều, đối chiếu mức xámbằng công cụ phần mềm để có thể thu nhận được ảnh chụp cắt lớp cuối cùng. Có thể nói đây là công đoạn mô phỏng được thực hiện thuận lợi nhất bởi vì Matlab là một công cụ phần mềm rất mạnh trong việc thực hiện các phép toán xử lý này. Mặt khác Matlab còn hỗ trợ rất nhiều các hàm xử lý về đồ hoạ cho phép hiển thị ngay được các dạng ảnh mong muốn. Theo công thức (1.37) tín hiệu suy giảm cảm ứng tự do (FID) thu nhận đươc sau bộ tách sóng pha cầu phương có dạng: ∫ +∆−−= dreerAtS rtriTt )]()([2)(.)( ϕωρ (3.16) Với |)(|* , rBA xyrOω= là hệ số tỷ lệ với độ nhạy của cuộn thu RF. Đối tượng )(rρ là mật độ proton tuỳ thuộc vào vị trí trong đối tượng. Từ các công thực (1.48), (1.49) về cặp biến đổi Fourier một chiều ta nhận thấy biểu thức (3.16) chính là biến đổi Fourier của hàm mật độ proton )(rρ với hàm trọng lượng là exp(-t/T2) và tham số của phép biến đổi là: ))()(( rtr ϕω +∆− . Như vậy thực hiện biến đổi ngược Fourier ta sẽ thu được ảnh mật độ proton )2/exp()( Ttr −ρ . Thực chất ma trận dữ liệu thô là một mảng dữ liệu hai chiều được thu nhận qua các bước mã hoá pha. Biểu diễn dưới dạng mô hình toán học của nó có dạng sau: 68 Với xGr Xγω =)( và PY yTGr γϕ =)( (3.18) Ở đây có 2 tham số thay đổi đó là tần số và pha được mã hoá theo vị trí x và y. Lúc này dùng phép biến đổi ngược Fourier 2 chiều ta sẽ thu được ảnh mật độ proton của cả một lát cắt. Trên cơ sở phân tích mô hình toán học ở trên ta tiến hành xây dựng sơ đồ thuật toán mô phỏng cho công đoạn này như sau: Hình 3.8: Sơ đồ thuật toán mô phỏng quá trình xử lý và hiển thị ảnh. Sau khi tải sẵn ma trận dữ liệu thô ta thực hiện xử lý bằng phép biến đổi Fourier 2 chiều để thu được ma trận ảnh. Khi đã nhận được ma trận ảnh, quá trình tiếp theo được thực hiện nhờ các hàm đồ hoạ để hiển thị được ảnh cuối cùng. Cũng có thể thực hiện tái tạo 3D với một tập hợp lát cắt nhất định hoặc thay đổi các tham số quan trọng của ảnh theo mong muốn. 69 3.2. Mô phỏng toàn bộ quy trình xử lý của thiết bị Trên cơ sở phân tích mô hình toán học và lưu đồ thuật toán mô phỏng của các công đoạn xử lý chính ở cả phần cứng lẫn phần mềm của thiết bị, ta tiến hành tổng hợp lại để đưa ra mô hình toán học và lưu đồ thuật toán mô phỏng cho toàn bộ thiết bị CCL CHTHN. Với việc mô phỏng toàn bộ quy trình xử lý, trước tiên ta thực hiện tạo giả một lát cắt đối tượng ở dạng một ma trận vuông (giả thiết ở đây là ma trận có kích thước NxN) Trong mỗi chu kỳ thu nhận, theo công thức (1.29) tín hiệu nhận được sẽ là tín hiệu tổng hợp của tất cả các tín hiệu phát ra từ mỗi voxel trong lát cắt: (3.19) )]().(2cos[].),(/exp[).,()( 2 ytxyxTtyxAtS Latcat ϕπω +−= ∑ Ở đây là biên độ của tín hiệu ứng với voxel tại toạ độ (x,y) trong lát cắt, giá trị của nó bằng giá trị của điểm (x,y) trong ma trận lát cắt nhập vào. Ở đây có thể hiểu chính là mật độ proton của voxel tại toạ độ (x,y) trong lát cắt tạo giả. ),( yxA ),( yxA Tương tự như vậy là hằng số thời gian dãn hồi của thành phần từ trường ngang ứng với tín hiệu phát ra từ voxel có toạ độ (x,y) trong lát cắt. ),(2 yxT )(xω là tần số của tín hiệu CHTHN tại các vị trí có toạ độ X. Đây là tần số đã được mã hoá bởi gradient mã hoá tần số, nó thay đổi tuyến tính theo chiều mã hoá tần số (chiều trục X) và được xác định theo công thức (1.47): xAx .)( EF0 += ωω (3.20) Ở đây AEF là hệ số mã hoá tần số, xGA .EF γ= (γ là hệ số từ quay, là gradient mã hoá tần số theo phương x) xG 0ω là tần số cộng hưởng gây ra bởi từ trường không đổi . 0B 70 )( yϕ là góc pha của tín hiệu CHTHN tại các vị trí có toạ độ Y. Đây là góc pha đã được mã hoá bởi gradient mã hoá pha, nó thay đổi tuyến tính theo chiều mã hoá pha (chiều trục Y) và được xác định theo công thức (1.48): yAy EP .)( 0 += ϕϕ (3.21) Ở đây AEP là hệ số mã hoá pha, τγ .. yEP GA = (γ là hệ số từ quay, là gradient mã hoá pha và yG τ là khoảng thời gian tác dụng của trường gradient này). 0ϕ là góc pha ban đầu của tín hiệu CHTHN khi thực hiện kích thích vô tuyến (với các phương pháp tạo ảnh có góc quay 900 thì 0ϕ =0). Điểm cần chú ý ở đây là thời gian t phải được lấy mẫu với tần số cao hơn ít nhất là gấp đôi tần số của tín hiệu CNHTHN để đảm bảo không xẩy ra hiện tượng nhiễu loạn hay “Aliasing”. Tín hiệu tổng hợp này được đưa đến bộ tách sóng pha cầu phương để loại bỏ thành phần tần số cao 0ω , đồng thời tách ra làm hai tín hiệu vuông pha nhau là phần thực và phần ảo như đã được biết ở các công thức (1.33) và (1.34). Phần thực: (3.22) )](...2cos[].),(/exp[).,()( EF2 ytxAyxTtyxAtS Latcat R ϕπ +−= ∑ Phần ảo: (3.23) )](...2sin[].),(/exp[).,()( EF2 ytxAyxTtyxAtS Latcat I ϕπ +−= ∑ Tổng hợp hai thành phần này lại dưới dạng phức: ))](...2(exp[].),(/exp[).,()( EF2 ytxAiyxTtyxAtS Latcat FID ϕπ +−−= ∑ (3.24) 71 Tín hiệu SFID(t) chính là tín hiệu suy giảm cảm ứng tự do (tín hiệu FID), tín hiệu này tạm thời được lưu trữ lại để tiếp tục thu nhận các tín hiệu ở các bước mã hoá pha tiếp theo. Để thu nhận được ma trận dữ liệu thô mang đầy đủ thông tin về vị trí và cường độ điểm ảnh ta cần phải thực hiện đủ số bước mã hoá pha cần thiết, cụ thể ở đây với ma trận lát cắt đối tượng có kích thước NxN thì số bước mã hoá pha chính bằng N. Khi đã thực hiện hết các bước mã hoá pha ta thu được ma trận dữ liệu thô với kích thước cũng bằng NxN. Thực hiện biến đổi ngược Fourier theo công thức (1.51) ta sẽ nhận được một ma trận ảnh với các điểm ảnh tương ứng với lát cắt đối tượng tạo giả ban đầu. Ma trận ảnh nhận được nhờ biến đổi ngược Fourier 2 chiều: ))](...2(exp[].),(/exp[).( )2( 1),( EF22 ytxAiyxTttSyxA Latcat ϕππ +−= ∑ (3.25) Bước tiếp theo là dùng các hàm đồ hoạ trong Matlab để xử lý và hiển thị ra ảnh cắt lớp mong muốn. Trên thực tế cần có một sự hiệu chỉnh nhỏ trên ma trận ảnh do có sự sai số, và dịch pha trong quá trình thực hiện phép biến đổi ngược Fourier 2 chiều. Dựa trên những phân tích về mô hình toán học ta tiến hành xây dựng sơ đồ thuật toán mô phỏng cho toàn bộ quy trình xử lý như sau: 72 Hình 3.9: Sơ đồ thuật toán mô phỏng toàn bộ quy trình xử lý Sơ đồ thuật toán mô phỏng này cho phép mô tả toàn bộ nguyên lý hoạt động của thiết bị CCL CHTHN theo đúng các quy trình mô phỏng ở từng công đoạn xử lý đã trình bày ở trên. 73 Tuy nhiên, như ta đã biết do có sự hạn chế về công cụ mô phỏng không cho phép thực hiện ở dải tần đúng như trên thực tế của thiết bị CCL CHTHN thật, nên ở đây khi thực hiện việc mô phỏng tổng hợp toàn bộ quy trình xử lý ta cần phải có một vài điều chỉnh cần thiết theo sơ đồ thuật toán sau: Hình 3.10: Sơ đồ thuật toán hiệu chỉnh lại 74 Điểm khác biệt trong sơ đồ thuật toán này so với sơ đồ trên là ở chỗ ta thực hiện mô phỏng bằng cách tách riêng vai trò của bộ tách sóng cầu phương. Khi này dạng công thức toán học của tín hiệu sẽ khác với dạng tín hiệu trong sơ đồ thuật toán trên, tức là có dạng phức như tín hiệu sau bộ tách sóng cầu phương theo công thức (3.24): (3.29) ))](...2(exp[].),(/exp[).,()( EF2 ytxAiyxTtyxAtS Latcat ϕπ +−−= ∑ Nguyên nhân của việc chỉ có thể thực hiện mô phỏng theo sơ đồ thuật toán này được giải thích cụ thể như sau: Như đã biết đối với thiết bị CCL CHTHN thật thì việc kích thích và thu nhận tín hiệu đều được thực hiện bằng phần cứng. Khi này tín hiệu CHTHN có tần số rất lớn ở dạng tương tự được đưa đến xử lý bởi bộ tách sóng cầu phương. Tín hiệu sau nó cũng là dạng tín hiệu tương tự nhưng có tần số thấp hơn rất nhiều, việc lấy mẫu tín hiệu này để chuyển sang dạng số không đòi hỏi tần số lấy mẫu phải quá cao. Khi đó có thể thực hiện lấy mẫu nó với số điểm tương đối thấp (khoảng vài trăm điểm) mà vẫn đảm bảo được tính chính xác của tín hiệu. Nhưng khi tất cả đều được mô phỏng trên máy tính, thì tín hiệu CHTHN cũng cần p