medical image analysis SVN (old)

#!/usr/bin/env python
#
# Copyright (c) Madrid 2008 
# BIT, ETSI Telecomunicacion, UPM
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program 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 General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#
#
#

from optparse import OptionParser

modules = {'Section' : [0, '', 'none://Section.py' ]
           }
import wx
import os
import sys
import string
from lxml import etree
from Section import *

def GetRelative(registered, segmented):
    if len(registered) != len(segmented):
        raise ValueError("both sets should have the same size")
    result = []
    
    for i in range(len(registered)):
        delta = registered[i] - segmented[i]
        #print " %f %f" % ( delta, segmented[i])
        result.append( ( 1.0 + delta * delta) / ( 1.0 + segmented[i] * segmented[i]) ) 

    return result


app = wx.App()


parser = OptionParser()

parser.add_option("-i", "--in-file-pattern", type="string", dest="inpat", 
                  action="store", help="input image file name pattern", default="reg%04d.exr")
parser.add_option("-s", "--segmentation-set", type="string", dest="segset", 
                  action="store", help="segmentation set file", default="segment.set")
parser.add_option("-r", "--reference-frame", type="int", dest="refframe", 
                  action="store", help="reference frame")
parser.add_option("-z", "--zoom", type="float", dest="zoom",
                  action="store", help="mask zoom", default=2.0)
parser.add_option("-t", "--shift", type="float", dest="shift",
                  action="store", help="intensity shift", default=20.0)


(options, args) = parser.parse_args()

file = open(options.segset, "r")
strset = file.read()
root = etree.XML(strset)
file.close()
set = WorkSet(node=root)

slice_nr = options.refframe

if slice_nr >=  set.GetFrameNumber():
    print "Frame number  %d not from range:[0-%d)" % (slice_nr, set.GetFrameNumber())
    exit(1)
            
masks = []
slice = set.GetFrame(slice_nr)
image = slice.GetImageZoomed(options.zoom)
size = image.GetSize()


for sec in slice.Sections:
    m = sec.GetMask(size, options.zoom)
    masks.append(m)

print "slice nr %d has %d sections" % (slice_nr, len(masks))

zw = size.x * options.zoom
zh = size.y * options.zoom

print "scale %d %d" % (zw, zh)
for i in range(set.GetFrameNumber()):
    frame = set.GetFrame(i)
    p_ref = frame.GetMaskedIntensities(options.zoom)

    image_file = options.inpat % (i)
    print image_file
    img = wx.Image(image_file, type=wx.BITMAP_TYPE_ANY, index = -1)
    if img.GetSize().x == 0:
        raise ValueError("Image file %s not found" % (image_file))

    img.Scale(zw,zh) 
    
    if len(p_ref) != 6:
        print "slice %d not segmented like it should be" % i
        exit(2) 
        
    p = []
    for mask in masks:
        val = 0.0
        n = len(mask)
        for m in mask:
            val = val + img.GetRed(m[0], m[1])
        p.append(val / n - options.shift)
    
    res = GetRelative(p, p_ref)

    print "%f %f %f %f %f %f" % (res[0], res[1], res[2], res[3], res[4], res[5])