지난번 올린 글에서 고민하던 voxelization algorithm이 해결되었다. 참으로 긴시간이었는데… 어찌보면 하이브리드 모델을 만드는데 있어서 가장 핵심적인 기술이라고 볼 수 있다. 이 부분에서 얼마나 정확하게 계산해내는가에 따라 하이브리드 모델의 장래가 결정된다.

처음 만들었던 voxelization algorithm의 핵심부분은 다음과 같다.
…..
        for t=0:step:1;
            x_pos=x(1)+(x(2)-x(1))*t;
            y_pos=y(1)+(y(2)-y(1))*t;
            voxel(xmin+round(x_pos/res(1)),ymin+round(y_pos/res(2)),zmin+round(z_slice/res(3)))=organtag(object);
        end

    temp=zeros(size(voxel,1),size(voxel,2),size(voxel,3)); % temporary matrix bounding object
    temp(find(voxel==organtag(object)))=1;
    for z=1:size(voxel,3)
    [B,temp(:,:,z)] = bwboundaries(temp(:,:,z),8,’holes’);
    end
    voxel(find(temp~=0))=0;
    voxel(find(imerode(temp,se)~=0))=organtag(object);
…..

이 코드의 문제는 본래 polygon mesh model의 volume을 너무 크게(20% 이상) 묘사한다는 것이다. ICRP89의 값에 3% 이내로 접근하려는 지금, 문제가 있는 코드였다. 이때부터 고민을 시작하여 다른 algorithm을 도입했다. 처음 것보다 속도는 훨씬 빨랐으나 여전히 volume문제가 있었다. 두번째 코드의 핵심부분이다.
…..
        % index voxels on the line segment
        if (x(2)-x(1))/(y(2)-y(1)) >= 0
            x_min=min(x); % slope>0
            x_max=max(x);
            sign=1;
        else
            x_min=max(x); % slope<0
            x_max=min(x);
            sign=-1;
        end
            y_min=min(y);
            y_max=max(y);
            if max(x)-min(x) > max(y)-min(y)
                % 0<slope<1
                for xtemp=x_min:sign*res(1)/2:x_max;
                    if y_max==y_min ytemp=y_min;
                    else
                    ytemp=((y_max-y_min)/(x_max-x_min))*(xtemp-x_min)+y_min;
                    end
                    voxel(xmargin+fix(xtemp/res(1)),ymargin+fix(ytemp/res(2)),zmargin+round(z_slice/res(3)))=organtag(object);
                end
            else
                % slope>1
                for ytemp=y_min:res(2):y_max;
                    if x_max==x_min xtemp=x_min;
                    else
                    xtemp=((x_max-x_min)/(y_max-y_min))*(ytemp-y_min)+x_min;
                    end
                    voxel(xmargin+fix(xtemp/res(1)),ymargin+fix(ytemp/res(2)),zmargin+round(z_slice/res(3)))=organtag(object);
                end
            end
…..

그냥 대충 봐도 뭔가 조잡한 느낌이다. 여전히 문제를 안고 있었고… 그래서 고민을 다시 시작했다. 결국에는 scanline algorithm이라는 것을 발견했고 오랜 코딩과 디버깅을 거쳐 결국에는 Voxelizer라는 코드를 완성했다. 이 코드의 핵심부분은 이곳에도 올려놓을 수 없다. 이전 코드의 문제점은 volume말고도 또 한가지가 있었는데 contour내부에 있는 contour는 묘사하지 못한다는 치명적인 오류였다. Cranium과 vertebrae쪽에서 문제가 많았다. Voxelizer는 이 모든 문제를 해결했다.

이 일은 내가 한 것이 아니라 하나님께서 하셨다는 것이 너무나 분명하다. 내 머리의 한계를 철저하게 경험한 뒤 기도하면서 한 일이라서 그렇다. 잠3:5-6말씀을 몇번이나 연습장에 쓰면서 했는지 모르겠다. 하나님께 모든 영광을 드리고 나중에라도 이 순간을 기억하며 교훈을 삼기위해 이곳에 글을 남긴다.