CAD CAM CNC Australia CAD Forum
CAD CAM Forum for support and discussion regarding AutoCAD, progeCAD IntelliCAD, Alibre, Mathcad, PTC Pro/E 2D 3D design T-FLEX and PartMaster CNC software
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Mathcad for Medical Imaging and Nuclear Medicine  ADD TO STUMBLEUPON

Post new topic   Reply to topic    CAD Forum Forum Index -> CAD and CAD/CAM Software
View previous topic :: View next topic  
Author Message

Joined: 12 Aug 2009
Posts: 6
Location: Australia

PostPosted: Mon Sep 14, 2009 5:54 am    Post subject: Mathcad for Medical Imaging and Nuclear Medicine Reply with quote

Here is another glimpse into our upcoming Mathcad webinar series. As usual, community and professional feedback is always welcome. We also welcome your suggestions about Mathcad applications in various industries. In Australia users can also test-drive our free Mathcad download for themselves.

Before we start, I also want to mention that we can be contacted in Australia with questions about Mathcad HERE.

With all radiologist and radiographers, one thing they all know that is important is that the quality of the images taken for the diagnostics must be of the highest quality possible. However, there are many factors that contribute to poor quality image acquisition such as parameters settings for the scan, image contrast, contrast sensitivity, distortion, noise as well as artifacts and blurring. Getting that balance between sensitivity and selectivity creates the need for image processing a necessity.

This Webinar will allow for image processing possible to be done in MathCAD, which can be used in the medical field where every little detail makes the difference in giving an accurate diagnosis. This webinar will be broken down into 5 parts as look at each process in MathCAD to give quality images without comprising the integrity of the scans taken.

1. Equalisation
2. Function and level mapping
3. Noise and Error measurement
4. Crisping
5. Filtering Noise

But why bother using Mathcad to do image processing when other Photoshop programs are available? The answer is simple. Mathcad allows for far greater control on how defined you would like the image to be without comprising the image as well as being able to customize using different algorithms. Also another factor to consider is that the images can be swapped in and out quickly to use the same function or that particular setting just by changing the address line where the image is located.

First of all, we will be using an image which has been provided through PTC for this tutorial. In the Handbook for image processing, the image which will be used to demonstrate the various features of the image processing pack in Mathcad is brain.gif as shown here. Through out this tutorial the original image is shown next to the edited image which has been enhanced by MathCAD with appropriate names.


Equalisation allows for the scanned image to be more defined by controlling on how the light and dark values are distributed in defined cumulative histogram of the image. This will in turn create a linear looking cumulative histogram of the scanned image and giving sharper details on the image.

Now to activate this function, command line typed out as equalize(M) , as M is defined as your image from the previous line of calculation to this command. However it is best if you use a similar layout to what I will be doing for the example of equalisation here:

Firstly, define the image in which is being read. The example show here is that I will be using M to define my image from the MathCAD image processing booklet.

M :=READ_IMAGE(“C:\Program Files\Mathcad\Mathcad 14\Handbook\improc\brain.gif”)

Next I will output the image to a histogram to see the degree of spread of the intensities in the 255 greyscale bands. Note that the number 256 is used since that the intensities start at and include 0 to 255, therefore giving 256 intensities). So the command of

H:=imhist(M ,256)

And for the histogram to work, we need to define the data in which the histogram will be using which is the pixel matrix of the image as we defined as so:

k:= 0..rows(H) – 1

Next create a Histogram and define the axis accordingly to the Hk as the function and k as the axis values for x and you should get a histogram like this:

Seeing the image is slightly dark, we will be spread out the intensities to define the features on the scan to help get a better detailed image of the scan. For this image, the cumulative histogram is given by the difference equation for C. So we will define as the following and repeat the same steps to create a cumulative histogram as shown with these steps.

Cumulative histogram:

Now that we can see the slope is not linear from the histogram, we will apply the equalise function here to see what happens. We do this by typing the following commands to generate the histogram and the new changes shown

Requal :=equalize(M)
J :=imhist(requal,256)

The new cumulative histogram of equalized image generated by typing the following commands:

Now just generate the cumulative Histogram:

Now as we can see the cumulative histogram is showing a relatively linear curve here give us the following images as a result.

Now as we can see the newly edited image has more defined features since that the equalisation has brought out the lines within the Cerebellum is more defined as well as the sulci and the gyri in the cerebral cortex.

CADDIT SMB & Enterprise CAD CAM Support
CAD CAM Software
Back to top
   View user's profile Send private message Visit poster's website

Joined: 12 Aug 2009
Posts: 6
Location: Australia

PostPosted: Mon Sep 21, 2009 4:45 am    Post subject: Function and Level Mapping for MathCAD Webinar Reply with quote

Function and Level Mapping

Function and level mapping shows the different levels of intensities across an existing image which will allow for different areas of the image to be more defined. This can be done to eliminate the amount of background noise as we saw in the equalisation step for the first webinar.

To activate this function, type in funmap(M,f) where M is the image matrix which needs to be generated first. The character f in the command is for the function to be performed at each vector or cells in the image matrix (256 levels to map the different intensities) which is going to be used to help define the details on the image.

This means that every time that the image matrix is processed by the new function, the mapping image matrix will be also calculated with the new changes to the image at each intensity each time separate as it is applied. To use this function, we type in

R := READ_IMAGE("C:\Program Files\Mathcad\Mathcad 14\Handbook\improc\brain.gif").

This allows us to setup first variable for the function mapping. Next we will use a function to help us generate the desired effect. Now there are a number of other functions which can be used for this or custom made for their desired effect. However this will require doing some experimenting and testing of the function applied to the image. For now we will be using the following function to create the desired effect.

Here is a small list of possible functions which can be used for the function mapping.

Once we have completed that we enter in

fmap := Re(funmap(R,F))

Since that we have generated the function as well the image matrix to use, using an output image you can see the results.

Level Mapping

Level Mapping allows for the replacement of the intensities with in a specified image by a specified area or vector of intensity. In other to simply put it, to increase the intensity levels within a specific area by use a defined vector. An example of this can be said to be the same of having a 29th element in a vector will give a new level for the pixels with an intensity of 29. It is important to note that images have entries within 0 and the length of the vector used of minus 1.

An example of this is that we would like to create an image with a squared intensity scale. We would create a Vector within:

r := 0 …255

This will result in the creation of the following curve constructed given us the what the vector will look like.

Now given by imaging pack we can use a number of examples as shown here, to refine the image for better screening and printing resolution.

This therefore, helps in enhancing the image for better diagnosis of the patient. Note to remember is that a particular function map or look-up table can be created with monitor or sensor, which maps irregularities cause by the display to their corrected values.

Once we are satisfied with the vector created we apply the vector to the Level mapping function with the following command and specify the image to be used in the level mapping as well.

R := READ_IMAGE("C:\Program Files\Mathcad\Mathcad 14\Handbook\improc\brain.gif").

level :=levelmap(R, vec)

Last edited by Senseineo on Mon Sep 21, 2009 6:06 am; edited 1 time in total
Back to top
   View user's profile Send private message Visit poster's website

Joined: 12 Aug 2009
Posts: 6
Location: Australia

PostPosted: Mon Sep 21, 2009 4:47 am    Post subject: Error and Noise Measurement for Mathcad Webinar Reply with quote

Error and Noise Measurement

With Error and Noise Measurement, we use functions which are based on the relative error( squared error ration, the mean squared error and the signal-to-noise ration between the two images which are used to be compared. These functions are used to determine the level of noise that affects an image after that it is processed or transmitted. To demonstrate this function we will be looking at a few examples here.

Firstly we need to define out variables here: which are R that represents out first or control image matrix that we are using. And Q which is the second image matrix , the same size as the first. Note that the functions return a number which represents the relative error, the mean squared error, or the signal-to-noise ratio (SNR) between M and Q. Remember that all returned values are in decibels (dB).

Now the first of the three commands that we are going to be looking at with error and noise measurement is the relative error. This function returns the squared error ratio over all the elements of the two matrices that are defined by M and Q. This is activated by the command of relerror(M, Q).

Now to workout the squared error ratio is as follows:

Over two matrices, this function rearranged and defined as:

After getting the required values for N by working out the ratios, simply type in the following command to work out the relative error for R:

err :=relerror(M ,N)
err :=0.0666667

To demonstrate the effect of this function, we will be applying this function to our example image here and define it:

R := READ_IMAGE("C:\Program Files\Mathcad\Mathcad 14\Handbook\improc\brain.gif").

Next we will add some random noise to this image by the following command and then calculate the relative error from the image matrix.

Q := addnoise(R,.05,128)
err :=relerror(R ,Q)
err :=0.059

Another way in which you can measure the amount of error in a image is the mean square error function. This gives the average squared error between the selected images of R and Q. For example we can read and define an image first. Once we have done that we can apply the re-quantize to 13 levels on it:

R := READ_IMAGE("C:\Program Files\Mathcad\Mathcad 14\Handbook\improc\brain.gif"
Q := imquant(R,13)

The mean squared error (MSE) is defined as follows:
I := 143
J := 234

Next we define the Dimensions of images:
I :=0..(I-1)
J :=0..(J-1)

MSE= 41.189

Calculating the MSE with the built-in function we obtain the same result:
MSE2 :=immse(R,Q)
MSE2 :=41.189

The third of the error and noise measurement functions in this section, is the SNR or the Signal to Noise Ratio. This command calculated the singal to noise ration of an image in decibels (dB). For example, let’s say when transmitting a scanned image between diagnosticians via a communication channel. The receive we observe the same image, but is corrupted by some noise interference. We will mimic or recreate this situation by applying the first scan as follows:

Q :=addnoise(R,0.3, 150)

The SNR is defined as the ration of the averages of power between the original and of the noise itself. The noise is obtained by subtracting the original image matrix with the recreated noisy image. To activate this function by typing in the function: SNR :=imsnr(R,Q)

This gives us a value of 3.57 (dB)

Now another commonly used function in image integrity is peak signal to noise ratio (PSNR). This function can be activated by using the following command:

Where A is defined by the image quality used (8 – 24bit, A = 255 for 8bit and A=4095 for 12 bit)
Back to top
   View user's profile Send private message Visit poster's website

Joined: 12 Aug 2009
Posts: 6
Location: Australia

PostPosted: Mon Sep 21, 2009 4:48 am    Post subject: Crisping for Mathcad Webinar Reply with quote


Sharpening the details of image is important to give definition to the image when diagnosing a patient. Mathcad provides crisping functions which help in this matter. They are as follows:

orthogonal crisping (R)
dia – crisping(R)
uni –crisping (R)
orthogonal5 crisping(R)

All of these functions work by the convolution of the specific crisping kernel within an selected matrix defined as R. The Number in the orthongal5 shows that the 5 x 5 kernel is used instead of a 3 x 3 kernel. Crisping can be used to restore lost sharpness to an image which as degraded due to transmission or image processing. The command will result in giving out a more defined or crisper image matrix. The edges remain unchanged since the kernel doesn’t overlap completely in these areas.

These are the matrix kernels that are used for the following functions

Now to use any of these functions, firstly define the image to be used and then define the command to be used. For Example, we will use the same scan through this webinar

R := READ_IMAGE("C:\Program Files\Mathcad\Mathcad 14\Handbook\improc\brain.gif"

Next we blur the image and crisp the image again.
R :=orthosmooth(R)
Orthogonal := orthocrisp(R)

Due to the nature of Crisping functions allowing full floating point calculations and compression of image values, we would need to rescale and equalize to spread the values out again. So in order to save time and space in the worksheet, we will combine these functions all into single line commands as shown here:

Orthogonal :=equalize(scale(orthogonal(R),0,255))

dia :=equalize(scale(diacrisp(R),0,255))

uni :=equalize(scale(unicrisp(R),0,255))

Ortho5 :=equalize(scale(orth5(R),0,255))

Here are some examples of combining commands and the results of crisping the original Image.

If you want to try Mathcad yourself, you can download the Mathcad trial HERE.

Last edited by Senseineo on Mon Sep 21, 2009 6:20 am; edited 2 times in total
Back to top
   View user's profile Send private message Visit poster's website

Joined: 12 Aug 2009
Posts: 6
Location: Australia

PostPosted: Mon Sep 21, 2009 6:07 am    Post subject: Noise Filtering for Mathcad Webinar Reply with quote

Noise Filtering

Noise filtering requires experience in the types of noise that the person is dealing with in each image. Depending on what generated the noise, what sort of edge distortions are acceptable in the output and what sort of noise it is will great depend on what sort of filtering which will be used. We will start of with a simple scattering of noise (otherwise known as salt and paper noise). We will simulate this and see if this can be reversed.

Orig :=READ_IMAGE(“C:\Program Files\Mathcad\Mathcad 14\Handbook\improc\brain.gif”)
Noisy: = addnoise(orig,.2,128)

Now to remove this type of noise there are two ways in which can be done to this. Firstly, the use of an averaging function such as smoothing or secondly a more complex method of median filtering.
So we will look at the first method of doing this the smoothing filtering. Smoothing of the image can be done by the following command:

Smooth := unismooth(noisy)

To use the median filtering, type in the following command:
Med :=medfilt(noisy)

Now the following images will show the difference of the filtering methods as we show the orginal image, the noise generated one, smoothing and the median filtering images.

As the images prove, the median filtering is a better filtering method for the scatter noise because it doesn’t change the intensity levels of the images, but uses what is available in the image. To observe these levels, 4 histograms are provide from each image to see the differences between them. To enable this simply type in

Intensity := imhist(<name of the image>,256)
k:= 0.. length(intensity1) – 1

As shown through the histogram, the original and the median filtering histogram exhibit very similar intensities proving that median filtering is a much more appropriate method in filtering the noise.

Median filtering is a better way to filter this type of noise because it does not change the intensity levels in the image, but only uses available levels. To see this, look at the histograms of the four images.

The next webinar will cover the wavelet smoothing for noise filtering.

If you want to try Mathcad yourself, you can download the Mathcad trial HERE.
Back to top
   View user's profile Send private message Visit poster's website
Display posts from previous:   
Post new topic   Reply to topic    CAD Forum Forum Index -> CAD and CAD/CAM Software All times are GMT
Page 1 of 1

Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum

Powered by phpBB 2.0.21-6 (Debian) © 2001, 2005 phpBB Group