00001 % DeOblique
00002 % a script that figures out the rotation to apply to
00003 % cardinal-plane-acquired data to bring it into alignment with oblique functional data
00004 %
00005 % Usage:
00006 % The script will prompt you for an Ifile filename, say: may0603mj/I.610
00007 % and will output the 3drotate command needed to apply to the
00008 % High resolution anatomical volume to bring it into alignment with
00009 % the Ifile (may0603mj/I.610)
00010 %
00011 % Hit enter to stop the loop. The results are logged in a file called
00012 % DeOblique.txt
00013 %
00014 %This script, with its silly name, was never meant to be distributed,
00015 %it was hastily written and meant for my own consumption only.
00016 %It has been only tested for slices slightly off the axial plane.
00017 %
00018 %If you are using this script for some reason, use it with
00019 %caution. Verify your results with AFNI at the end.
00020 %
00021 %ziad@nih.gov
00022
00023 FuncName = 'DeOblique';
00024
00025 %load image header
00026 %RAS here means LPI
00027
00028 %fname = sprintf('/home/ziad/tmp/pm040803/5to10/023/I.138');
00029 %fname = sprintf('/home/ziad/tmp/pm040803/6to10/024/I.138');
00030 %fname = sprintf('/home/ziad/tmp/pm040803/spgr/005/I.073');
00031
00032
00033
00034 fname = input ('Enter I filename: ', 's');
00035 while (~isempty(fname)),
00036 [su_hdr,ex_hdr,se_hdr,im_hdr,pix_hdr,im_offset] = GE_readHeader(fname);
00037
00038 %get three points forming corners of plane and the center of the plane
00039 tlhc1 = [im_hdr.tlhc_R im_hdr.tlhc_A im_hdr.tlhc_S];
00040 trhc1 = [im_hdr.trhc_R im_hdr.trhc_A im_hdr.trhc_S];
00041 brhc1 = [im_hdr.brhc_R im_hdr.brhc_A im_hdr.brhc_S];
00042 ctr1 = [im_hdr.ctr_R im_hdr.ctr_A im_hdr.ctr_S];
00043 Tri1.XYZ = [tlhc1; trhc1; brhc1];
00044 [err,Eq1] = Plane_Equation (Tri1);
00045
00046 %an idea of the displacement at the edges, typically the minimum of d
00047 d = 2 .* (tlhc1 - ctr1);
00048
00049 %form the equation of the desired plane
00050 ctr2 = ctr1;
00051
00052 CardPlane = sprintf('Axial');
00053 %CardPlane = sprintf('Sagittal');
00054 %CardPlane = sprintf('Coronal');
00055
00056 %desired plane is one that is parallel to Axial plane and passing through ctr2
00057 switch (CardPlane(1)),
00058 case 'A',
00059 % Axial plane eqation is z - cnst = 0;
00060 Eq2 = [0 0 1 -ctr2(3)];
00061 case 'C',
00062 % Coronal plane equation is y - cnst = 0;
00063 Eq2 = [0 1 0 -ctr2(3)];
00064 case 'S',
00065 % Coronal plane equation is x - cnst = 0;
00066 Eq2 = [1 0 0 -ctr2(3)];
00067 otherwise,
00068 fprintf (1,'Error Bad Plane.\n');
00069 end
00070
00071 % find the intersection line of the two planes
00072 u = cross(Eq1(1:3),Eq2(1:3)); %the direction vector of the intersection line
00073 %should check that cross product is not ~ 0 (i.e.
00074 u = u ./ norm(u, 2);
00075
00076 %This line must pass by ctr because I said so
00077 somescale = abs(max(trhc1 - brhc1))./3;
00078 Line = [ctr1; ctr1+somescale.*u];
00079
00080
00081 %Good, now get the angle
00082 CosAngle = dot(Eq1(1:3),Eq2(1:3)) ./ ( norm(Eq1(1:3), 2) .* norm(Eq2(1:3), 2));
00083 AngleRad = acos(CosAngle)
00084 AngleDeg = AngleRad .* 180 ./ pi
00085
00086 %Show the results
00087 Opt.Fig = figure(1); clf
00088
00089 subplot 211; cla
00090 plot3(Tri1.XYZ(:,1), Tri1.XYZ(:, 2), Tri1.XYZ(:,3),'r*'); axis tight;hold on
00091 plot3(ctr1(1), ctr1(2) , ctr1(3), 'b*');
00092
00093 [err,PatchHandles] = ShowPlane (Eq1, Opt); view(3); hold on
00094 [err,PatchHandles] = ShowPlane (Eq2, Opt); view(3); hold on
00095 xlabel ('R'); ylabel ('A'); zlabel ('I');
00096
00097 plot3(Line(:,1), Line(:,2), Line(:,3), 'r-', 'LineWidth', 4);
00098 stmp = sprintf('Rotation about intersection line from %s plane to oblique is %g degrees.\nMax displacement from cardinal plane %g mm\n', ...
00099 CardPlane, AngleDeg, min(d));
00100
00101 title (stmp);
00102 subplot 212; cla
00103 plot3(Tri1.XYZ(:,1), Tri1.XYZ(:, 2), Tri1.XYZ(:,3),'r*'); axis tight;hold on
00104 plot3(ctr1(1), ctr1(2) , ctr1(3), 'b*');
00105
00106 [err,PatchHandles] = ShowPlane (Eq1, Opt); view(3); hold on
00107 [err,PatchHandles] = ShowPlane (Eq2, Opt); view(3); hold on
00108 xlabel ('Left to Right'); ylabel ('Posterior to Anterior'); zlabel ('Inferior to Superior');
00109
00110 plot3(Line(:,1), Line(:,2), Line(:,3), 'r-', 'LineWidth', 4);
00111 axis equal; hold on
00112
00113 %suggest a rotation with 3dvolreg, but first, make sure this is not a double oblique slice
00114 N_RotAx = 0
00115 if (Line(1,1) ~= Line (2,1)) ,
00116 RotAx(1) = 1;
00117 N_RotAx = N_RotAx + 1;
00118 else
00119 RotAx(1) = 0;
00120 end
00121 if (Line(1,2) ~= Line (2,2)) ,
00122 RotAx(2) = 1;
00123 N_RotAx = N_RotAx + 1;
00124 else
00125 RotAx(2) = 0;
00126 end
00127 if (Line(1,3) ~= Line (2,3)) ,
00128 RotAx(3) = 1;
00129 N_RotAx = N_RotAx + 1;
00130 else
00131 RotAx(3) = 0;
00132 end
00133
00134 if (N_RotAx > 1),
00135 fprintf (1,'Looks like you have super oblique slices, not ready to deal with that yet.\n');
00136 else
00137 if (RotAx(1)),
00138 AlignToOblique = sprintf ('3drotate -rotate %gR 0A 0I', AngleDeg);
00139 end
00140 if (RotAx(2)),
00141 AlignToOblique = sprintf ('3drotate -rotate 0R %gA 0I', AngleDeg);
00142 end
00143 if (RotAx(3)),
00144 AlignToOblique = sprintf ('3drotate -rotate 0R 0A %gI ', AngleDeg);
00145 end
00146
00147 fidv(1) = 1;
00148 fidv(2) = fopen('DeOblique_out.txt','a');
00149 for (outp=1:1:2),
00150 %suggest 3drotate command
00151 fprintf (fidv(outp),'\nTo rotate Data to match %s:\n%s -prefix Data_Rotated_to_Oblique Data\nHit Enter to continue\n',...
00152 fname, AlignToOblique);
00153 %if (outp == 1) pause; end
00154 fprintf (fidv(outp),'\nTypically, Data represents an anatomical \n');
00155 fprintf (fidv(outp),' data set acquired in a cardinal plane \n');
00156 fprintf (fidv(outp),' (Axial, Sagittal or Coronal). The oblique slice\n');
00157 fprintf (fidv(outp),' is from functional data. The suggested 3drotate\n');
00158 fprintf (fidv(outp),' command is meant to bring the anatomy (Data) into\n');
00159 fprintf (fidv(outp),' alignment with the function.\n');
00160 fprintf (fidv(outp),'If the rotation appears to be in the wrong direction,\n');
00161 fprintf (fidv(outp),' try using %g for a rotation angle.\n', -AngleDeg);
00162 fprintf (fidv(outp),' When that occurs, please notify me. \n\tziad@nih.gov\n\n');
00163 end
00164 end
00165
00166 if (0), %other untested methods ....
00167 % ----------------------------------
00168 % get the rotation matrix
00169 %Making Obliques Cardinals:
00170 [err,ObliqueNowCard.XYZ, rot1] = AxisRotate3D(Tri1.XYZ, u , AngleRad, ctr1);
00171 for (i=1:1:3),
00172 fprintf (1,'Oblique p%g: | %.2f %.2f %.2f |\t%s p%g: | %.2f %.2f %.2f |\n', ...
00173 i, Tri1.XYZ(i,1), Tri1.XYZ(i,2), Tri1.XYZ(i,3), ...
00174 CardPlane, i, ObliqueNowCard.XYZ(i,1), ObliqueNowCard.XYZ(i,2), ObliqueNowCard.XYZ(i,3));
00175 end
00176 %Tri1.XYZ
00177 %ObliqueNowCard.XYZ
00178 %rot1
00179
00180 %Making Cardinals Obliques
00181 [err,CardNowOblique.XYZ, rot2] = AxisRotate3D(ObliqueNowCard.XYZ, u , -AngleRad, ctr1);
00182 for (i=1:1:3),
00183 fprintf (1,'Cardinal p%g: | %.2f %.2f %.2f |\tObique p%g: | %.2f %.2f %.2f |\n', ...
00184 i, ObliqueNowCard.XYZ(i,1), ObliqueNowCard.XYZ(i,2), ObliqueNowCard.XYZ(i,3), ...
00185 i, CardNowOblique.XYZ(i,1), CardNowOblique.XYZ(i,2), CardNowOblique.XYZ(i,3));
00186 end
00187
00188 %ObliqueNowCard.XYZ
00189 %CardNowOblique.XYZ
00190 %rot2
00191
00192 %write results for 3drotate.
00193 % ************* SWAP rot1 and rot2 to make rtation go in the correct direction ...
00194 fmatvec = sprintf ('%s_To%s.matvec', fname, CardPlane);
00195 fid = fopen (fmatvec, 'w');
00196 if (fid < 0),
00197 fprintf (1,'Error: Failed to open output file %s for writing.\nCheck your permissions.\n', fmatvec);
00198 end
00199 for (i=1:1:3),
00200 fprintf (fid, '%g %g %g 0\n', rot2(i,1), rot2(i,2), rot2(i,3));
00201 end
00202 fclose (fid);
00203 fprintf (1,'To turn the oblique data to %s plane use:\n', CardPlane);
00204 fprintf (1,'3drotate -matvec_dicom %s ...\n\n', fmatvec );
00205
00206
00207 [err, PathString, FileString] = GetPath(fname);
00208 fmatvec = sprintf ('%sCardinal_To%s.matvec', PathString, FileString);
00209 fid = fopen (fmatvec, 'w');
00210 if (fid < 0),
00211 fprintf (1,'Error: Failed to open output file %s for writing.\nCheck your permissions.\n', fmatvec);
00212 end
00213 for (i=1:1:3),
00214 fprintf (fid, '%g %g %g 0\n', rot1(i,1), rot1(i,2), rot1(i,3));
00215 end
00216 fclose (fid);
00217 fprintf (1,'To turn the cardinal data to %s''s plane use:\n', fname);
00218 fprintf (1,'3drotate -matvec_dicom %s ...\n\n', fmatvec );
00219
00220 fprintf (1,'To test for individual points add options: -origin %g %g %g -points\n\n',...
00221 ctr1(1), ctr1(2), ctr1(3));
00222
00223 end
00224
00225 fname = input ('Enter I filename: ', 's');
00226 end
00227 fprintf (fidv(2), '\n\n\tCreated with %s\n\t%s, %s\n', FuncName, pwd, date);
00228 fclose(fidv(2));
